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ABSTRACT 


Acoustic ray paths in the ocean are known to exhibit sig- 
nificant horizontal deflections after repeated reflections 
from the bottom. The effect may be quantitatively and quali- 

. tatively observed through a ray trace model which permits a 
change in direction of the vertical plane of propagation as a 
function of bottom slope and grazing angle. ECTRACE is a 
family of computer programs which traces a bundle of rays in 
three dimensions and utilizes bottom depth values as a portion 
of its input. Included are related programs which develop sea 
bed models from digitized bathymetry data or synthetic bathymetry 


functions. 
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I. INTRODUCTION 


ECTRACE is a FORTRAN computer program designed for the 
investigation of bathymetric effects on horizontal acoustic 
ray tracing. Its supporting computer programs include routines 
for constructing a grid model of the sea bed from contour data, 
ray tracing routines, and associated printed and plotted output. 
The programs were written for IBM 360/370 series systems with 
a Versatec plotter and related software. The significant fea- 
ture of the tracing algorithm is its ability to shift the direc- 
tion of the vertical plane of propagation after a bottom 
reflection. A sea bed modeled by triangular facets fitted 
through three points of a grid cell makes this feature possible 
with a minimum effect on computation time. 

Several models for tracing acoustic rays in the ocean are 
used in the Navy today. These models are based on a simplified 
environment because of practical limitations on data and com- 
puter run time, and they adequately trace a bundle of rays in 
a single vertical plane. Some enable intensity calculations 
at a given distance from the source. MTRIMAIN [1], for example, 
is a model which features range as well as depth dependence 
of sound speed and accounts for some interaction with irregu- 
lar bathymetry. These models have a large number of applica- 
tions in deep ocean cases involving ducts, channels, and 


convergence zones, and bottom reflections where the plane of 


the ray path does not exhibit horizontal deflection. 


Actual cases of long range propagation over irregular bot- 
toms show any single ray path rarely remains within a single 
vertical plane. Any application of ray theory to propagation 
involving bottom interaction must account for a horizontal 
deflection caused by ray heading changes after repeated reflec- 
tion from a tilted bottom plane. This effect is most pro- 
nounced for long range paths over smoothly sloping surfaces, 
but any valid ray approximation can be seen to undergo signif- 
cant heading changes after only a few reflections from an 
irregular or undulating bottom. 

The behavior of sound propagation over troughs, ridges, 
wedges, and seamounts has been studied analytically by Har- 
rison [2] using normal mode theory and ray invariants. Since 
actual ocean bottom topography over large area defies prac- 
tical analytical description, except in the stochastic sense, 
DeWitt [3] developed numerical techniques for three-dimensional 
ray tracing in which the bottom consists of triangular facets 
generated numerically from actual or interpolated bathymetry 
data. This approximation scheme greatly simplified the prob- 
lem of calculating the three-dimensional ray path parameters 
after a bottom reflection and provides a means of constructing 
reasonably accurate ray traces (for a verification, see Appen- 
dix F). 

The ECTRACE modeling program is an adaptation of the 
Dewitt technique for use on an IBM 360/370 series computer. 
The initial programs may be used to generate synthetic bot- 


tom topographical grids or create approximations to real 


ocean bottom grids from a bathymetric data source file. The 


topographical features of the generated grid can be depicted 
by a contour plot and a perspective surface plot for compari- 
son with the original contours. The primary program is the 
ray tracing program which accepts the generated bottom grid 
and performs stepwise ray position calculations from user- 
selected initial heading and elevation angles. The output 

is the printed history of each ray, and a two-dimensional plot 
of the ray paths projected on a vertical plane intersecting 
the bottom profile along the axis of the ray fan. Addition- 
ally, a horizontal plane projection is plotted which may be 
overlaid on the appropriate contour plot to correlate the 
horizontal curvature effect with grid bathymetry. Punched 
card output from several ray tracing runs may be combined in 
one independent routine for a comprehensive visual inspection 
of acoustic shadow zones in the horizontal plane. 

It is very important to note that ECTRACE permits refrac- 
tion caused by sound velocity gradients in the vertical plane 
only and that the sound velocity in the horizontal plane is 
assumed constant within a well defined water mass. All hori- 
zontal bending effects revealed by ECTRACE are the result of 
bottom reflection only. 

With the ability to convert actual bathymetry data to a 
matrix of depth values, discrete faceted approximations of a 
selected oceanic region can be used for generating ray traces 
of operational significance. In addition, qualitative studies 
of horizontal ray deflections over idealized bottom configura- 
tions are enhanced by the computer-generated plots of ray paths 


over easily modifiable synthetic bathymetry. ECTRACE and the 
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augmenting programs are envisioned to serve as a basis for 
future development of three-dimensional ray modeling and inves- 


tigation efforts. 


15 


II. DESCRIPTION OF ECTRACE AND PERIPHERAL PROGRAMS 


This chapter describes each of the programs, their pri-~ 
mary subroutines, and their interrelationship in producing a 
ray trace which may be used to investigate the bathymetric 
effects on horizontal ray curvature. They are presented in 
their intended order of use. Flow Diagram (Fig. 1) and Table 


I are provided to assist in understanding program relationships. 


A. GENBOT 

This program converts contour data into a matrix of depth 
values for storage on a permanent device. It also makes a 
two dimensional contour plot of the input data and super- 
imposes reference latitude and longitude lines (Fig. 2). 

The contour plot can be placed underneath the horizontal ray 
trace plot of ECTRACE or ECCOM on a light table to aid in 
the visual investigation of the ray curvature (see Figs. 3 
and 4, and Figs. 5, 6, and 7). 

The bathymetry contour data used in this research exists 
as a sequence of data sets representing ten by ten degree 
regions of the North Atlantic Ocean and were obtained exclu- 
sively from a tape compiled in a joint project of the Naval 
Research Laboratory and the Mobil Oil Company. In their 
original form the data were sufficient only for producing a 
computer plot of region contours, thereby reproducing por- 


tions of a chart titled "Bathymetry of the Norwegian-Green- 


land and Western Barents Sea" [{4] (see Fig. 8). Their use 


in GENBOT required us to insert the depth values of the con- 
tour lines and transfer the results to a separate tape (de- 
scribed in Appendix A). 

The following subrovtines are called by GENBOT: 

1) CORNER calculates boundary latitudes and longitudes for 
screening only those data needed for interpolation. 
GEOPLT superimposes reference latitude and longitude 
lines on the subregion contours. 

GEODST calculates the forward or inverse solutions of 
the geodetic triangle to accomplish the above. 

RAIN 1, RAIN 2, and RAIN 3 [5] were obtained from the 
SSP3 program library at NPS and are used to perform 
interpolations on data points for construction of the 


depth matrix. 


B. SYNGEN 


This program produces a synthetic bathymetry grid for 


storage on a permanent device. It may be used to generate 
the following types of ideal sea bed configurations: wedge, 
trough, ridge, conical seamount (Fig. 9), sinusoidal undula- 
tion, or parabolic basin (Fig. 10). These grids can be used 
in the ECTRACE program as a depth matrix in the same manner 
as the grids of actual bathymetry produced by GENBOT. SYNGEN 


is discussed in detail in Appendix B. 


C. G3DP 
Program G3PD produces a perspective of any portion of the 
generated depth matrix through use of the NPS system subroutine 


CONTUR ({4]. Figures 9, 10, and 11 are examples. G3DP uses 
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job control commands to link with the GRDSCT subroutine of the 


ECTRACE load module to extract the desired area (hereafter 
called the working matrix) to be plotted. Appendix C describes 


the application of G3DP in detail. 


D. GRDCHK 

This program is used to check the validity of the depth 
matrix generated by the programs GENBOT and SYNGEN. It produces 
a contour plot and a vertical plot of user-selected rows of the 
generated matrix. By comparing consecutive rows the existence 
and location of anomalous data points may be determined (Figs. 
4 and 12). It contains the previously mentioned subroutine 
GEODST, CORNER, and a version of GEOPLT (called GEOPLR) to draw 
reference geographic lines for comparison of the two contour 
plots (GENBOT vs GRDCHK). Refer to Appendix D for additional 


discussion of GRDCHK, 


E. ECTRACE 

The ECTRACE program is a package of subroutines which 
exist in load module form on a permanent storage device in 
the NPS computer system user library. Its use requires a 
user supplied calling program which calls subroutine TRACER 
and contains the necessary job control language (JCL) to link 
with the module, These procedures are explained in more de- 
tail in Chapter IV of this report. 

ECTRACE traces rays in three dimensions. The sound speed 
field is assumed to vary piecewise-linearly with depth only, 


yet provisions are made to permit simulation of up the three 


distinct water masses separated discontinuously by vertical 
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fronts. The bathymetry structure, generated by one of the 
two grid-making programs, is a three dimensional surface 
modeled by discrete triangular facets fitted through cells 
of adjacent depth values. Rays are traced as piecewise arc 
segments each with a radius of curvature dependent on the 
vertically-structured sound speed gradient. 

When the rays undergo a bottom bounce, specular reflection 
is assumed. Bottom loss values can be calculated from the 
Subprogram FBLOSS (described later), a standard Navy model 
such as the FACT bottom-loss routine, or a user supplied 
routine, The trace history supplied by ECTRACE includes the 
location and characteristics of each ray reversal (surface 
reflections, bottom bounces, and refractive turning points). 
The printed data lists the new ray parameters determined 
after the reversal had occurred. Chapter V explains the 
printed output in more detail. 

The ECTRACE plot product includes a plane view of the 
ray paths projected on a vertical plane, a plot of sound 
speed profiles, and a horizontal (plane) view of the ray 
paths (Figs. 5 and 13). The norizontal view can be overlaid 
on the appropriate contour plot generated by GENBOT or GRDCHK, 
as in Figs. 12 and 14. 

The following subroutines are used by ECTRACE: 

1) Tracer is the primary tracing subroutine which calls 
all other subroutines of the ECTRACE program and is 
capable of processing a bundle of rays separated into 


ray fans. A ray fan is defined here as a number of 
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rays (user determined) each with a different initial 


heading measured clockwise from grid north, but all 


with the same initial elevation angle. A positive 


elevation angle describes a ray which is pointed toward 


increasing depth. 


GRDSCT reads into core a fixed-dimension section of a 


depth matrix file of arbitrary size. The ray trace is 


restricted to the boundaries of the core-loaded working 


matrix. A ray which departs the working matrix boundaries 


is terminated and program control is passed to the next 


ray's trace. 


IDTSUB identifies the three matrix subscript pairs whose 


depth values define a plane triangular bottom facet 


underneath the ray segment head. 


DEEP calculates the depth at any horizontal position by 


solving for the z-coordinate on a plane fitted through 


the three depth matrix values. 


CONTAC calculates the three-dimensional coordinates of 


the intersection of the ray path and triangular bottom 


facet, 


BOANG calculates the grazing angle, elevation angle, and 


new heading of a bottom reflected ray. 


FBLOSS returns bottom loss values in dB as a function 


of grazing angle for a bottom reflected ray. These 


values are based on NRL geophysical survey data of the 


Greenland-Norwegian Sea and Iceland environment collected 


in the mid-1970's and have been widely averaged for 


convenience, The user has the option of substituting a 
more specialized function if bottom loss values are 
critical. 

IDPROF determines the water mass, and hence the sound 
speed profile, affecting the ray with each segment 
iteration, 

NUPROF initializes calculations for new water mass param- 
eters once the ray has passed the boundary. 

CHNLIM determines in advance the refractive sound channel 
limits for each ray from its invariant. This is used to 
test a ray at its turn-around point and to reveal its 
type (surface reflected, refracted-surface reflected, 
bottom reflected-surface reflected) in the printed 

trace history. 

ANGPRT is used to print a summary of selected bottom 
contact in parameters for each ray fan. 

The following subroutines each have only one calling 
statement in TRACER. The subroutines they in turn invoke 
are dependent upon the computer graphic system installed 
and would require internal modification if exported to 

a system without Versatec software. For additional 


information NPS users should consult Ref. 7. 


a) BGNPLT initializes the plotter and draws all borders, 


titles, and axes. 
BDTPLT draws a vertical profile of the sea bed cen- 
tered along the mean heading of the ray bundle, as 


shown in the upper plot of Figure 13. 
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RNGPLT traces all ray segment projections on the 


vertical plot, 

T2DPLT plots the horizontal track of each ray. Each 
ray fan uses one of thirteen symbols to indicate the 
bottom bounce positions of each ray. Figure 5 shows 
some of the available symbols. 

SSPPLT draws the sound speed profiles. 

ENDPLT draws the horizontal plot legend and termi- 


nates all plotting. 


F. ECCOM 

Program ECCOM uses the optionally produced punched card 
output of one or more ECTRACE jobs to draw a composite hori- 
zontal plot of several ray bundles, enabling a comprehensive 
pictorial study of shadow zones and acoustic convergence in 
the horizontal plane. The composite horizontal plot is de- 
signed to portray the rays which emanate from the same source 
or from multiple sources along a desired track. The ray num- 
bers assigned and punched for each ray on each ECTRACE run 
assist the user in singling out rays of interest (or non- 
interest) before combination. Figures 7 and 18 are examples 
of this combining technique. 

Amplifying remarks on these programs and subroutines can 
be found in the comments in the computer FORTRAN listings. 
For additional description of plots and computer printed out- 


puts, refer to the appropriate appendices. 


III. PROGRAM LOGIC AND THEORY 


ECTRACE traces the three-dimensional paths of a bundle of 
rays, one ray at a time. The bundle is divided into fans of 
rays of the same initial elevation angle. A separate ray fan 
is traced for each fixed increment of initial elevation angle 
between the limits ELST and ELEND (input constants). A posi- 
tive elevation angle describes a ray pointed downward. The 
initial ray headings in each ran fan extend from HDST to HDEND 
(input constants), measured in degrees clockwise from grid 


north. 


A. BOTTOM AND SOUND FIELD MODELING 
Figure III.1 shows the projection of the bottom surface 


onto the horizontal x-y plane. 
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FIGURE III.1 Bottom Surface Projection 


The projection of each triangular facet is a right isosce- 
les triangle. The entire surface can be defined by specifying 
the depth values at the vertices of the projection triangl?2s. 
ECTRACE stores these depth values (km) in the two-dimensional 
working array ZB. The parameter DHN (km) represents the spac- 
ing between the vertices in both the X (East) and Y (North) 
directions. Thus, the value of the matrix point ZB(I,J) is 
the depth at the vertex point with X coordinate (I - 1) * 

DHN and Y coordinate (J - 1) * DHN, relative to the ZB origin. 

The matrix values for ZB are read from a source file stored 
on a permanent storage device such as a disk pack. A matrix of 
dimensions 151 by 151 with spacing DHN = km has been used for 
trial runs. 


1. Bottom Contact Point 


The subroutine CONTAC locates the coordinates and calcu- 


lates the depth of the bottom bounce point. The parameters 
describing locations of the ray segment head and tail are among 
those passed to the main program. In this subroutine, the ray 
segment is treated as unrefracted (a straight line), a reasonable 
approximation for the small distances involved in determining 

the contact point. 

First the horizontal positions of the ray segment head (x,y) 
called CEE and CNN in CONTAC, and (x,,y,) and called CRE and 
CRN in CONTAC, are checked to see if they both lie in the same 
projection triangle. If not, an interative routine searches 
for two points on the segment such that the horizontal coordi- 


nates of both points lie inside the same triangular cell and 


define the portion of the segment that penetrates the bottom 
facet. Then the depth of the ray head z, and the depth of the 
tail z, are computed. 

Figure III.2. illustrates a profile of a cell bottom with 
a ray intersecting it in the propagation plane. The solution 
of the triad (Kor Q0Z,) identifying the point of bottom contact 
proceeds as follows: 
Using subroutine DEEP (described in Section 3) compute HS and 


H 


27 the bottom depth at (x, -¥,) and (X,Y), respectively. Then 


compute: 


where (XQ 1YQ12Q) are contact point parameters, 


z= H,Z, - HZ, 


co be 
Az + H, - H, 
If az #0,M= 77 4, 
Az 
If az =0,M= 71 7 Hy 
H, - H, 
Then xX, =x + Mr sin 6 and 
Yo =y +Mr cos ¢ 


where > is the ray heading (angle between the y-z plane and the 
propagation plane). 
2. Projection Triangle Vertex Subscripts 
The subroutine IDTSUB calculates the I subscripts, IRT, 


ILT, IVR and the J subscripts JBT, JTP, JVR for the vertices of 
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FIGURE III.2. Profile of the bottom cell in the propagation 
plane, with a ray contact, 

of the projection triangle which contains the ray's horizontal 

coordinates x and y. The subscript parameters are indicated 


in Figure III.3. The following suffixes identify the relative 


locations of the vertices of the projection triangle: 


RT - right 
LT - left 
TP - top 


Lat, ITP 
BT - bottom 


run, Ter IRT, ToT 
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FIGURE III.3. Projection Triangle Vertices 


3. Depth Calculation 
Subroutine DEEP calculates the depth in kilometers 


corresponding to the horizontal coordinates x and y. Auto- 
matically, IDTSUB is called to identify the projection tri- 
angle subscripts. The depth Za is calculated by solving the 
equation of a plane defined by the depths at the vertex points 
(D., D,, and D.,) Figure III.4 illustrates the calculations 
made based on the orientation of the projection triangle out- 
lining the bottom facet plane. The quantities Xo and Y, repre- 
sent the horizontal position (in units of DHN) of the depth 
calculation point relative to the horizontal position of the 
depth calculation point relative to the horizontal position 

of the depth value D, (the lower left vertex of the projec- 
tion triangle) as calculated in Secion 1. The depth values 

Ds. and DY are normalized by the length of the equilateral 


side 1, in units of DHN. Finally, the depth is solved by 


Zo = xD. + Y.>, + D, 
B. RAY LOCATION AND DIRECTION PARAMETERS 

A three dimensional left hand coordinate system is used 
(Fig. III.5.), with the positive z axis pointing in the direc- 
tion of increasing depth and the sea surface level lying in 
the x-y plane. As the ray is traced through three dimensional 
space, the program variables which are used to locate the tail 
of the ray vector are: 

DEP1 - depth coordinate 


CRE - x coordinate in km 
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FIGURE III.4. Depth calculation on bottom facet 
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CRN - y coordinate in km 
CRE and CRN are components of the horizontal range increment 
DR which forms the propagation plane heading angle PHI with 


the y-axis. 
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FIGURE III.5. Ray Location Parameters 


Given that the ray vector originates at the point with 
the above coordinates, the direction of the ray vector is 
specified by the following angles: 

TH] - ray elevation angle in the plane of propagation 
measured from the ray vector's projection on the xy 

plane to the ray vector. A ray pointing toward the 

ocean bottom will have a positive elevation angle. 

Theta takes on values between -1/2 and +7/2 (output 
values are in degrees). 

PHI ~ ray (propagation plane) heading angle between the 
positive y-axis and the projection of the ray vector on 
the xy plane, measured clockwise from the positive y-axis 
(grid North). PHI takes on values between -7 and f. 


Initial values for TH1 and PHI are calculated from user input 


values (in degrees), which are converted to radians for compu- 


tational purposes. The unit vectors in the horizontal plane 


are given as UVX = sin (PHI) and UVY = cos (PHI). 


C. RAY PATH CALCULATIONS 
1. Path Length 

ECTRACE confines refraction to the vertical plane alone. 
It is recognized that while horizontal sound speed gradients 
exist, their effects are sufficiently small within a well defined 
water mass to be neglected compared to the bathymetric effect. 

However, while the horizontal gradients are likely to 
be slight, important discontinuities in water characteristics 


may occur horizontally across oceanic fronts. For this reason, 


ECTRACE has been designed to simulate up to three distinct 
water masses and two associated frontal boundaries. 

Within a water mass having uniform characteristics 
along fixed depths a ray path conveniently defines a vertical 


plane within which Snell's Law can be written as 


V = c(z) /cosé (z) (1) 


where V is treated as invariant within the water mass (termed 
the vertex velocity in some texts since it is the speed at the 
depth at which the ray vector becomes horizontal). The term 
c(z) is the sound speed at depth z and 6(z) is the elevation 
angle at depth z. 


Equation (1) uses the simplification that the speed of 


propagation and the elevation angle are functions of depth only. 


It is generally possible to define a stratified medium in which 
each layer's gradient is constant. Although such a model re- 
quires a large number of data points to adequately describe a 
water column of complicated sound speed structure, the piece- 
wise linear gradient approximation greatly simplifies the ray 
path descriptions. 

Taking advantage of the approximations, we use a simple 


equation for the sound speed gradient, within a layer, 


g = dc/dz (2) 


Integrating yields 


where the constant of integration is avoided by arbitrarily 
specifying the origin of the coordinate axes at the depth where 
c = 0. By applying Snell's Law we obtain 

Vv 

Zz =—— cos 

5 (4) 

and 
dz --= sin6d6 (5) 


Specifying r as the horizontal distance axis of the vertical 


plane and applying analytical geometry there results 


dz = tanédr (6) 
dr = -— cos6dé (7) 
and 
V : 
r=- — siné 
g (8) 


when measuring r from the vertical (z-axis) to the point at 
which 6= 0 (horizontal projection). 


Squaring (4) and (8) and summing yields 


yr? - 2? = (—~) (9) 


This is the equation of a circle of radius (V/g) and whose 
center is at the origin of the coordinate system just described. 


The equation 


eee 
er (10) 


describes the radius of curvature of a ray path in a constant 
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sound speed gradient. The minus sign has been chosen to allow 
p to be positive when the refracted ray is increasing its ele- 
vation angle. 

The ccordinate system and computer-approximated ray 


path segments are shown in Fig. III.6. 


FIGURE III.6. Ray Path Geometry 


It can be shown from Fig. III.6. that the ray segment 


path length between P and Q is 

AL = pA8s (11) 
where 

A6 = 6, = 8, (12) 


Since the initial elevation angle of the ray is presumed known, 


successive values of @ are found iteratively: 


REPS rere Tr “Teer 


Simple geometry also allows for solution of the depth change, 
Az = = ole ta p(cost - cos6 ) ' (14) 

and likewise the increase in horizontal range, 
Ar zr - 5. p(sinde | - siné ) : (15) 


Note that in the Fig. III.6., A&, Az, and Ar must all be positive, 
so the sign of o must agree with that of A686. Equation (10), 
which gives o the opposite sign of the gradient, satisfied this 
requirement. Only Az is allowed to become negative, as when a 
ray decreases its depth. 
2. Travel Time 
From the diagram (Fig. III.7.) the sound speed relation- 


ship becomes 


c= (16) 
dz 
dt = Gsing (17) 


Zr = on: ea 


z x 
FIGURE III.7. Ray path and sound speed relationship. 
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The ray invariant 


(18) 


and the trigonometric identify sin6 = \ - cos’6 combine 


and give dt as a function of depth 


dt = nS nee (19) 


e(t) f = fac(z)1)* 


The time-of-flight of the acoustic ray along the path is found 


by integrating from the initial point to the final point: 


daz 


f pee eee 
i C(2) in = factz)) au 


Before integrating it is assumed that the gradient is constant 
along the segment. More generally, ECTRACE assumes a constant 
gradient in the layer bounded by Zz and zoe 


Then 
c(z) = c(z) +t g(z - z ) 
1 2 i 
where z sz s zo 
The integral may now be evaluated by introducing anew variable 
W= such that 
W=2 
Then dW = dz and C(W) = gW. Therefore, 


aw 


gw fp - a2g*w? 


Ee ee 
eye ey (23) 
Wek +01 = a-g-WE-] 
W. (1 + sind.) 
a i Zn [ bt See earn (24) 
g W, (1 + sind,) : 


Values for W are restricted to being in the same layer and on 


the same side of a turning point in the ray path. 


D. RAY TRACING PROCEDURE 

ECTRACE begins tracing a ray from its given initial coordi- 
nates CRE, CRN, and DEP1l. First the water mass is identified 
to determine the local gradient and therefore the ray invariant 
and radius of curvature. The ray path is constructed in the 
form of arc segments using an increment DTH calculated from a 
chosen arc length DL. Then the new 6 value is calculated and 
Equations (14) and (15) are used to determine the ray segment's 
depth increment and horizontal range increment. The horizontal 
position is updated using the horizontal unit direction vectors 
(UVX and UVY) calculated from PHI, the initial heading of the 
ray's plane of propagation. 

The interim update 3-D coordinates of the ray head may now 


be stated in FORTRAN as 


CEE = CRE + DR * UVX (25) 
CNN = CRN + DR * UVY (26) 
DEP2 = DEP1 and DZ (27) 


where DR represents the horizontal range increment Ar and D2 is 


the depth incrementdAz. 
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E. RAY CONSTRAINT CONDITIONS 

Each time the ray path is incremented, certain constraint 
conditions must be checked. These constraints are bottom re- 
flections and the sound channel limits. 

1. Bottom Reflections 

After the interim horizontal coordinates CEE and CNN 

are calculated, they are passed to DEEP to find the bottom 
depth WH2 at the ray head. If DEP2 is greater than or equal to 
WH2, a bottom bounce has occurred. In this case the values of 
DEP2, CEE, and CNN are adjusted by linear interpolation and a 
point of bottom contact is established. Subroutine BOANG is 
then called to calculate new values of TH1, PHI, UVY, and the 


ray grazing angle GRAZ. 


A 
a 


bottom facet 


Plane of 
propagation 
+ intersection 


FIGURE III.8. Bottom Bounce Reflection Angles. 


The basic calculations of BOANG involve the vectors: 


=> 


- unit ray vector before bounce 


> 


- unit vector normal to triangle facet at the bounce point 
vi- unit ray vector after bounce 

The components of v are calculated from TH1 and PHI. The com- 

ponents of h are calculated from the equation of the plane of 


the bottom facet. The new vector ¥ is then calculated from 
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the vector equation, 
visv-2(ven)n, (28) 


New values of TH1, PHI, UVX and UVY are then calculated from 
the compnents of vi. Finally, the grazing angle GRAZ comes 


from 
sin(GRAZ) =-(v +n) . (29) 


Reflection from a sloping bottom plane facet, as op- 
posed to a surface reflection, forces a reevaluation of the 
Snell's Law constant V, and therefore the radius of curvature i 
R. The invariant V remains in effect until the next bottom 
bounce occurs or a new water mass is entered. Thus, after 
summoning BOANG, the ray tracing continues with the new ray 
parameters. The sound channel limits CHS and CHD described 
below are also recalculated at this point since they also 
depend on V, and GRAZ is passed to function FBLOSS to calcu- 
late a bottom loss value. 

2. Sound Channel Limits 

The Snell's Law equation (1) says that the ray trace 
is constrained to depths where the speed of sound does not 
exceed V. A given ray path may be constrained in the sense 
that it is purely refracted (R), refracted or reflected (RSR 
or RBR) or purely reflected (BRSR). The determining factors 
are V and the local sound speed by profile. V is calculated 
from the sound speed at the initial depth and the initial 


elevation angle, and is recalculated after a bottom bounce or 


upon entering a new water mass. 
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CHS and CHD are the refractive sound channel limits 
at the shallow and deep extents respectively. Since the sound 
speed profile is approximated by piecewise continuous linear 
gradients, the water column is treated as a column of indexed 
layers, each characterized by upper and lower (shallow and 
deep) depth boundaries and a constant sound speed gradient. 
The limiting depths CHS and CHD are found by scanning the 
water mass layers both above and below the current ray loca-~ 
tion until a boundary sound speed value is found which exceeds 


V. The limiting depth is then calculated from 


Zam = 2 + (V -c)/g (30) 


where z and c are evaluated at the top of the layer and g is 
the gradient of the limiting layer. 

For cases where g = 0 or where V is not exceeded by 
any sound speed in the local profile, CHS and CHD are set to 


the physical limits of the water column. 


F. RAY PROPAGATION TERMINATION CONDITIONS 
Tests are made during each ray segment iteration for events 
which terminate the ray trace. These events are: 
* Escape from grid - when the ray departs the finite x and 
y constraints of the working depth array (ZB). 
« Depth boundaries - when the depth detected at the ray is 
less than ten wavelengths or greater than 10 km, 


* Bottom loss - when the attenuation due to bottom reflec- 


tion exceeds a specified input amount in dB. 


* Error at bottom contact - when defective program code 
or anomalous depth data prevent the interpolative rou- 
tine CONTAC from recovering the bottom depth within a 
reasonable number (50) of iterative attempts. 

* Trapping - when it becomes apparent that the ray will 
not contact tire bottom, rendering further three-dimen- 
sional tracing unnecessary. 

* Bottom bounce gate - when the total number of bottom 
reflections exceeds a vser specified amount. 

The cause of ray termination is specified in the printed 

ray history, giving the ray parameters at the time of 


termination, 


G. IMPROVEMENTS MADE OVER THE ORIGINAL NRL VERSION 

The ray tracing algorithm originally existed as a program 
called ABOUNC at the Naval Research Laboratory. ECTRACE dif- 
fers significantly from ABOUNC in the following areas: 

1. Ray Path Development 

ABOUNC calculates ray parameters through right tri- 

angle solutions and fnell's Law, but develops the ray path 
along fixed depth. Since this prevents rays from becoming 
horizontal, a reversal is forced when the magnitude of the 
elevation angle goes below a gate value. In contrast, ECTRACE 
develops the ray path along are segments and radii of curva- 
ture, fixing only the maximum length of the segments while 
interrupting the increment as necessary when a segment reaches 
a reversal point or layer boundary, or to keep the increment 


within a maximum amount. 
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2. Computer Word Storage 
ABOUNC floating point variables are nearly all FORTRAN 


double precision, The largest single burden on computer core t 


serene et 


is the depth matrix, read from a magnetic tape in card image 
format, each card record containing eight items (fields) ten 
bytes long. This format permitted a depth range of +9999 km 
with resolution of +l cm. By reducing the depth range to +999 
km and resolution to +10 cm, two unnecessary bytes are removed 
from each field, enabling a fit of ten items per record and a 
per.aanent storage reduction of twenty percent. In addition, 
ECTRACE reads the matrix into a single precision array to re- 
duce the core requirements for the data by fifty percent. 
3. Grid File Manipulation 

ABOUNC reads the contents of the entire depth matrix 
data file into array whose dimensions must be exactly equal to 
those of the matrix. This requirement is not inconvenient on 
a large computer with FORTRAN dialect which permits objects 
time dimensioning. ECTRACE is written in FORTRAN G for an IBM 
Series/360 Model 67 machine whose limitations required some 
fundamental modifications of the program and its job control. 

To maintain flexibility and reduce run time, ECTRACE 
was reduced to a package of subprograms which have been pre- 
compiled and stored on a permanent device in load module form. 


The user must supply a calling program which serves primarily 


to establish the dimensions of the depth array. The user's 
selection of array dimensions are made on the basis of core 


economy and area of interest rather than the dimensions of the 
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input matrix. ECTRACE, through its internal subroutine GRDSCT, 
extracts a rectangular working subregion from any part of the 
input matrix. Trial runs of ECTRACE used a working array of 
dimensions 151 by 151 representing a sector of 150 km on a side 
with a core requirement of less than 250K, while the input ma- 
trix was 369 by 443. Without GRDSCT, ECTRACE would have re- 
quired an allocation of 800K bytes using the same input file. 
4, Treatment of Small Elevation Angles 

As stated previously, ABOUNC forced a reversal of the 
ray as the elevation angle of the leading segment approached 
a gate value close to zero. Since the ECTRACE tracing scheme 
permits a full range of elevation angles, early trial runs 
revealed an inherent error phenomenon as @ approached zero in 
a steep gradient. Specifically, a ag on the order of 3 de- 
grees caused by a small radius of curvature or along arc seg- 
ment traced a ray to its vertex depth in error as much as 50 
meters from the Snell's Law prediction. Since the depth change 
was calculated from the difference in elevation angle cosines, 
the error was found to arise from the machine inability to 
retain difference precision for any pair of numbers close to 
integer values (round off error). Since the difference rather 


than the actual cosine values is needed for this calculation, 


2 p% 
the small angle series approximation cos 1 - 5 + 5 - a ue 8 
was used to formulate the equation 
cosé - cosé = 5(67~@7) = 4(6% -6"), (31) 
1 2 ; ar a4 201 


When restricted to small angles, this approximation yields 
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greater accuracy than the straightforward machine calculation 
of the left hand side, since the precision retained by the 
registers is very high. 
5. Travel Time 
ABOUNC calculates the travel time of an acoustic sig- 
nal along a ray segment by dividing the linearly-approximated 
path length by the mean sound speed. ECTRACE accounts for 
path curvature and depth dependence of the sound speed through 
integral evaluation, discussed in detail in section C.2. 
6. Track Events 
Layer transitions, reflections, and turnarounds are 
events which may require a change in ray parameters. ABOUNC 
tests for these events in algorithmic sequence, makes the 
necessary changes in response to the first event detected 
and continues the trace from the current position of the ray 
vector head. ECTRACE retraces a ray segment to the event 
depth and tests for all other events before beginning the next 
segment from the transition point. This logic eliminates the 
possibility of a long ray segment crossing more than one event 


depth with only one event detected. 
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IV. USER INSTRUCTIONS (ECTRACE) 


A single step ECTRACE job traces a bundle of rays from a 
single source in the form of ray fans. Each ray fan is a bun- 
dle of rays of common initial elevation angle. All ray fans 
are bounded between to initial heading values measured from 
grid North (grid North equals true North at the grid center). 

During the following discussion the reader may wish to 
refer to Fortran listing in Appendix M. The user instructions 
for the other programs (GENBOT, G3DP, SYNGEN, GRDCHK, and 
ECCOM) are explained in their individual appendices (User In- 


structions and Output Description). 


A. JOB CONTROL CARDS (JCL) 

The following JCL cards and parameters are required for 

an ECTRACE run: 

* JOB card - The CPU time parameter should allow 60 sec- 
conds for every seven rays to be traced. 

« EXEC card - This card must specify the FORTCLGW proce- 
dure (Appendix 0). 250K bytes of core is required for 
the standard 22801 (151 by 151) element depth matrix 
portion plus 4K for every additional 1000 elements. 

The computer systems job class definitions should be 
considered before deciding upon the size of the working 
matrix and the number of rays. For example, 

EXEC FORTCLGW, REGION.GO=250K 


is a class C job at NPS when limited to five minutes of 
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time, which is adequate for tracing 40 rays. 
Calling program - Following the 
//FORT.SYSIN DD * 

card, the calling program must contain cards 

COMMON/DIMS/KMX ,KMY and 

CALL TRACER (ZB) . 
ZB is a REAL * 4 array dimensioned exactly by KMX,KMY, and 
these values (KMX,KMY) are initialized to the values 
of the DIMENSION statement using assignment statements, 
as in the following example. 

DIMENSION 2B(151,151) 

COMMON /DIMS/ KMX,KMY 

KMX=151 


KMY=151 


CALL TRACER (ZB) 


STOP 

END 
KMX and KMY are INTEGER * 4 and should be less than or 
equal to the dimensions of the data matrix from which 
the depth values will be taken. 
LINK cards - These cards point to the ECTRACE load 
module. They must include all members of the module 
for which the user does not supply a substitute. An 
example which includes all members is: 

//LINK.USDD UNIT=3330, VOL=SER=DISKO1,DSN=1270.ECTRACE, 

// DISP=SHR,LABEL=(,,,IN) 


//ULINK.SYSIN DD * 


INCLUDE USDD (TRACER, GRDSCT , IDPROF ,NUPROF , CHNLIM) 
INCLUDE USDD(BGNPLT ,SSPPLT, BDTPLT , RNGPLT , T2DPLT, ENDPLT) 
INCLUDE USDD(IDTSUB , DEEP, CONTAC , BOANG , ANGPRT , FBLOSS) 
ENTRY MAIN | 
/* 


¢ GO cards ~- ECTRACE uses data set reference number one 


for input of the depth matrix. For example, to point 

to a depth matrix residing in a data set called $9999.GRID 

on DISKO03, the following DD statement is used: 
//GO.FTOLF001 DD UNIT=330, VOL=SER=DISK03,DSN=S9999, 

GRID, DISP=SHR, LABEL=(,,,1IN) 

If a depth matrix from an outside source is used, it 

must be in the format specified in the Output Grid Data 

Format section of Appendix A. If a punched card output 

is not desired to make an ECCOM run, include the card 
//GO.FPT0O7F001 DD DUMMY ° 


ahead of the GO.FTOILFO01 DD statement shown above. 


B. DATA CARDS 
Following the 
//GO.SYSIN DD * 
card, the user selected data must be provided in the following 
format. FORTRAN field descriptors are in parentheses. The 
FORTRAN variable names are listed for reference, 
¢ Card 1 (20A4) 
OTIT - title of the run 
¢ Card 2 (958.2) 


XO - source position for x coordinate (kilometers). 
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source position y coordinate (kilometers). 
source depth (meters). A source depth which 
exceeds the water depth will be decreased to one 
meter above the bottom. Source coordinates ina 
nonpropagating location will also be adjusted. 
receiver position x coordinate (kilometers). 
receiver position y coordinate (kilometers). 
receiver depth (meters). The receiver may be 
positioned at the source if calcuations are not 
desired. Excessive depth is automatically 
corrected. 
frequency (Hz). A ray is terminated when the 
water depth becomes less than ten wavelengths 
during propagation. 
(optional) segment length (meters). The parameter 
affects the smoothness of the vertical plot. 
Smaller values increase smoothness and run time. 
Large values create the risk of missing a botton 
peak. The default value is one wavelength. 

Card (6F8.2) 

XMIN initial x value of horizontal plot (kilometers). 

YMIN - initial y value of horizontal plot (kilometers). 

The above values are used for drawing the axis of the 

horizontal plot and for determining the section of the 

input matrix to be read into a working depth matrix (ZB). 


The values for KMX and KMY, recorded cell resolution 


(DHN), and the limits of the input matrix ultimately 
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determine the axis ranges. The user should be Fanivias 
with these parameters when choosing the origin of the 
plot, 

zST - initial depth value for vertical plot (kilometers) 


ZEND - final depth value for vertical plot (kilometers) 
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ARST - initial range value for vertical plot (kilometers) 
AREND- final range value for vertical plot (kilometers) 
The above values are used to draw the axes of the vertical 

plot. The initial values may be negative. The program 
adjusts the final values if necessary for plot esthetics. 
* Card 4 (6F8.2) 
HDST - heading of the first ray in each ray fan (degrees) 


HDEND heading of the last ray in each ray fan (degrees) 


ELST - elevation angle of the first ray fan (degrees) 

ELEND - elevation angle of the last ray fan (degrees) 
Headings are measured from grid north. Elevations are 
measured from the horizontal plane, positive downward 
(toward increasing depth). 

HORLIM - maximum ducting or channeling range (kilometers). 
This is a gate value, like BLMAX CARD 2, which 
terminates the trace of a ray when exceeded. 
HORLIM serves to screen out those rays which 
do not appear to contact the bottom. 

SMNGL ~- (optional) small angle value for alternative 


4 cosine difference calculation (degrees). A 
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value of less than three degrees is recommended. 


See the discussion on "small angle approximations" 
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in Chapter III. The default value of two 


degrees is set when left blank. 


Card 5 (615) 
NHD - number of rays per ray fan (21) 
NEL - number of ray fans per bundle(2l) 


NSTOP - maximum number of bottom bounces per ray (21) 


The program restricts NHD to 30 rays, with a maximum 
total number of 100 rays produced, each with a maximum 


number of 50 bottom bounces or 200 reversals. 

INCPLT - number of ray segments defining one plot seg- 
ment (21). With sufficiently small ray seg- 
ments it is not necessary to use each iteration 
for plot definition to achieve a smooth plot, 
so this parameter is provided to enable the 
user to coordinate FDL and INCPLT for run 
time economy. 

NLEN - length of the longest horizontal plot axis 
(<19 inches). This value controls the physi- 
cal size of the plot. A value less than four 
will supress the plot (resulting in an ABEND 
code in the HASP log), but it will not affect 
the calculations or printed and punched out- 
put. A value greater than 15 requires addi- 
tional JCL for strip plotting, since the 
total plot width is NLEN plus six inches. 


Refer to Ref. 7 for strip plotting. 
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NSSP - number of sound speed profiles (water masses 
provided as data (1 < NSSP <¢ 3). 

Card 6 (A4,I3) 

OMAS (1) - name of first water mass 

NVAL (1) - number of (depth, sound speed) pairs defining 

the sound speed profile (1 < NVAL < 50). 

Card 7 - (through card 16 if necessary) (10F8.2) 
These are the data pairs (depth, sound speed), in meters 
an m/sec., from the sea surface to the last data point 
in the profile. The gradients will be calculated linearly 
between these points, with the final gradient below the 
last point set to 0.01668 m/sec. by the program. 
Next card (4F8,2) 
If the previous data is to be followed by the data for 
another water mass, this card must be included to define 
the frontal boundary between the two. The order in 
which the water masses appear as inptut determines their 
positions relative to the fronts. Cards which describe 
the next water mass, follow the frontal definition card 
and are formatted as described cards 6 and beyond. The 
last data card is the last string of (depth, sound speed) 
pairs of the last water mass; refer to Appendix M (sec- 
ond page) for an example of three water masses. The 
following variables define the boundary between each 
water mass; 

FENTX (1) - x coordinate of front beginning point (km) 


FENTY (1) - y coordinate of front beginning point (km) 
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FENTX (1) - x coordinate of front end point (km) 

FEXTY (1) - y coordinate of front end point (km) 
The examples shown in Fig. IV.1 identify the water 
mass boundaries for various input axis values in the 
form (FENTX, FENTY), (FENTX, FEXTY) for each water mass 
boundary. The first example a. represents the data 
from an original ECTRACE run, refer to the second page 
of Appendix M. Exambles b. and d. show the priority 
of the water masses. The first water mass defined, 
Northwestern Mohns Ridge (NWMR), dominates the second, 
Southeastern Mohns Ridge (SEMR), and the third, isove- 


locity (ISOV), where they overlap. 
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(0,46) (150,46) 
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FIGURE IV.1. Examples of water mass boundaries for various values of 
(FENTX,FENTY) and (FEXTX,FEXTY). 
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V. INTERPRETATION OF ECTRACE OUTPUT 


The following discussion refers to an ECTRACE run titled 
"Trial 67***Synthetic Basin***" (Pigure 5), which is a trace 
conducted in a synthetic parabolic basin (Figure 6). A por- 
tion of the printed output is shown in Tables II - V. In this 
ECTRACE run, eight ray fans were selected to be projected be- 
tween -30 and +30 degrees elevation. Each ray fan was com- 
posed of only one ray, thus the initial heading for each ray 


was 010 degrees. 


A. PRINTED OUTPUT 
Refer to Table II through v (Appendix G) as examples of 


the following discussion, 


1. Parameter Tableau (Table ITI) 


After processing and adjusting the input parameters 
just described, ECTRACE prints a final tableau of these parame- 
ters plus some additional parameters calculated from the origi- 
nal input such as layer gradients, the total number of rays 
to be traced, and the spacing between rays. In addition, the 
results of the depth matrix-grid load is summarized, revealing 
data from the header record of the grid file. 

2. Ray History Tableau (Table III) 

During the trace of each individual ray, this tableau 
prints a line of ray parameters at each end point and reversal. 
The following is a description of one line of parameters. 


* TYPE identifies the nature of the reversal, and is 


TRONS OT ATE Lea ORT TOS 


abbreviatec INIT for initial position, BOTT for bottom 
bounce, SURF for surface reflection, OVER for turnover, 
and STOP for terminated position. 

X, Y, and DEPTH are the coordinates in kilometers of the 
reversal and end point. 

GRID HEADING is the ray heading measured clockwise from 
grid North at the completion of the reversal. 

ELEV is the ray's elevation angle measured downward 
(positive toward increasing depth) from the horizontal 
zero degree reference). 

BRPNG TO REC is the bearing from the reversal point to 
the receiver position. 

DIST TO REC is the slant range distance from the reversal 
point to the receiver. 

WTR MASS is the water mass identifier indicating the 
sound speed profile applied in the present calculations. 
LYR is the constant gradient layer number corresponding 
to the order input as profile data (i.e., layer 1 is the 
surface layer and indicates that the surface gradient 
applies. 

DEPTH LIMITS are the physical and refractive sound chan- 
nel limits of the ray, determined by its invariant. If 
MIN is 0.0, an upward traveling ray will reflect off 

the surface before reaching its upper vertex (turnover) 
depth. If the MAX depth is greater than or equal to 10km, 
a line of asterisks is printed since the ray could never 


attain its turnunder depth before reaching the bottom. 
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Otherwise, these limits are the ray's vertex depths, 


and if the reversal type is OVER or UNDR, the appropri- 


re 


ate limit should exactly match the depth value on the 
same line. 

C is the sound speed at the reversal point. This is the ' 
ray invariant used for tracing. 

GRAZ ANGLE is the angle in degrees between the ray vector 
and the line of intersection of the bottom facet plane 
with the vertical plane. This is the angle used to 
calcualte bottom loss and the new ray direction vector. 
BTM ANGLE is the angle between the bottom facet plane 
normal and the vertical (z-axis), in degrees. 

BTM LOSS is the accumulated intensity loss in dB after 
the bottom reflection. When this value equals or exceeds 
BLMAX, the ray's trace will terminate. 

LPS is the number of steps required by the CONTAC sub- 
routine .7 locate the point of bottom contact, A long 
ray segment length may require more steps than a shorter 
one. A maximum number of 50 steps is permitted before 
the bottom facet is assumed anomalous and the ray's 

trace is terminated. 

Prior to the final line of an individual ray's trace 
history, is a summary of plot definition parameters plus 
data on the ray's CPA to the receiver, and reason for 


ray trace termination. The CPA calcualtions are accurate 


vA 


to within one ray segment length. 
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3. Ray Fan Summary Tableau (Table IV) 


At the completion of the trace of one complete ray 
fan, a table is printed listing each ray's elevation angle 
and depth in meters before each bottom bounce. All zero values 
at the bottom of each ray's column indicate blank data. 
4. CPA Summary Tableau (Table v) 
When all ray traces have been completed, a table is 
printed listing each ray's heading, time, and distance at CPA 


to the receiver. 


B. PLOTTED OUTPUT 
ECTRACE produces a single mc .tag? of up to five separate 


plots (Fig. 13, the total dimensions of which depend upon the 


é 


user's selection for the maximum size in incnes of the longest 


axis of the horizontal plot (NLEN). The following table lists 


some sample values for NLEN and the resulting dimensions of 
the compiete plot. 


Maximum Maximum 


NLEN (in.)Height Width Remarks 
5 7 8.33 minimum 
10 14 16.67 
15 21 25 maximum size without strip 
plotting [5] 
19 26.6 31.67 maximum allowable 


Components of the ECTRACE plot follow. 
1. Sound Speed Profiles (Figure 13) 
A plot of sound speed versus depth is made for each 


water mass defined. 
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2. Vertical Plot (Figs. 15 and 16) 


Most two-dimensional ray trace programs produce plots 
of ray paths on graphs of depth versus range. The vertical 
plot is similar except that all ray paths are projected from 
three dimensions onto the vertical plane. The plane is oriented 
along the mean ray fan heading and includes its profile of the 
sea bed, In ECTRACE trial 65B (Figure 13) the ray fans each 
contain two rays whose initial headings vary from 360 to 020 
degrees as indicated by 10.000 +/- 10.000 on the plot. The 
bathymetry profile is made along the heading of 010 degrees. 
On a plot, the heading value (HEADINGS) always indicates the 
direction of the bathymetry slice. The (+/-) value signifies 
the initial heading of the first and last rays and that addi- 
tional rays are evenly spaced in the interval. 

Since the rays are actually traced in three dimensions, 
bottom reflections occurring outside the vertical plot plane 
may appear to be taking place above or below the bottom sur- 
face instead of at the boundary (Fig. 15). Only those bottom 
contact points which occur in the projection plane must lie 
on the seabed profile. 

The symbols listed on the plot (Fig. 5) are used to 
identify bottom contact locations in the vertical and horizon- 
tal plot, and are unique for each ray fan except that up to 
thirteen symbols are used before repetition, as shown in 


Table VI. 


3. Horizontal Plot (Fig. 17) 
A 


This plot is the most important plot product of 
ECTRACE and is the plot of the ray paths on a horizontal plane 
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(topview). This plot displays the xX-y positions of each 
bottom contact point using the ray's (ray fan) identification 
symbol and connects them with straight lines, revealing the 


horizontal deflection effect. Each ray's plot begins at the 


source, ends up at the last bottom contact point in the grid 
and is labeled with the ray number (Fig. 17). The ray num- 
ber can be used to locate its printed and punched output 
history. 

The axes are labeled with respect to the origin of 
the input matrix. The geographic center position and coor- 
dinates of the input matrix, if included in the bounds of 
the working matrix (FORTRAN array ZB), is also plotted for 
reference, It is important to note that this plot, without 
the bottom contact symbols, may combine the results of other 
ECTRACE runs by combining the optional punched card outputs 


as input to ECCOM (Figs. 7 and 18). 


C. PUNCHED OUTPUT 
The punched card output for a single ECTRACE job consists 


of the following cards: 


CARD DATA CONTENTS 
1 The runtitle 
2-3 The grid title 
4 Grid reference coordinates and scale 
factor 
5 Initial elevation, initial heading, 


number of points (BB) and ray num- 
ber of the first ray 
6+ X-Y coordinate pairs of plot points 
of the first ray 
(Data for additional rays is repeated as 5 and 6+ above) 
Last 9999., to flag the end of data 
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These cards may be stacked as a member with those of 
several ECTRACE jobs as input to a single ECCOM job. For 
ECCOM to be meaningful, cards 2 through 4 of each stack mem- 
ber should be identical, and the user should verify this 
himself. 

Punched card output may be quite voluminous and should 
not be produced unless required. The following card, inserted 
before the first GO card (JCL), will suppress the punched 
output: 


//GO.FTO7F001 DD DUMMY . 
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4 VI. SUMMARY 


A. POTENTIAL APPLICATIONS 

ECTRACE, in its present form, has many potential 

applications: 

« The investigation of asympotic deflection angles in 
ocean areas where bathymetry data are well-documented, 
allows optimum hydrophone positioning. For example, 
Ref. 6 cites cases where seemingly attractive receiver 
positions would be affected by shadowing in the hori- 
zontal plane. 

« Adaptation of ECTRACE ray tracing methods could be 
incorporated into existing programs which presently 
rely upon an assumption of radial symmetry of ray 
propagation about a source or receiver (no horizontal 
deflection). 

« The horizontal deflection effect may lend additional 
insight into the study of travel times and intensities 
of acoustic signals, enabling more accurate source 


fixing information. 


‘ B. IMPROVEMENT 
There are many areas of improvement which would increase 


the utility of ECTRACE and its supporting computer programs 


for the above applications: 
* A better routine to convert bathymetry contour data 


into a matrix of depth values is essential for modeling 
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regions for which only contour data exists. Most known 
interpolation algorithms perform best on data which are 
randomly or regularly spaced rather than in the form of 
contours. For example RAIN2 [5] would probably perform 
better with the data used to make the contours than it 
does with the contours themselves. While the contour 
data are useful for cartographic and geologic purposes, 
bathymetry data in a less processed form are significantly 
better sources for constructing computer models of the 
sea bed. 

ECTRACE in its present form makes no calculations of 
intensities or their losses due to spreading, attenua- 
tion, or scattering, thus preventing its use as a 
propagation loss model, 

An improved bottom loss function is recommended over 
the supplied FORTRAN function saibprogram FBLOSS. The 
user may perform this substitution easily by deleting 
FBLOSS from the INCLUDE step in the LINK statement of 
the ECTRACE job control and supplying another function 
subprogram (added at the end of the calling program) 
that accepts a double-precision argument for grazing 
angle and is also called FBLOSS. 

The assumption of specular reflection neglects the 
phenomenon of transmission into the sediment and its 
effect on the location of the reflected ray due to sub- 
bottom boundaries. It would be interesting to incorporate 


in an improved version of ECTRACE. 
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APPENDIX A 
GENBOT USER INSTRUCTIONS AND OUTPUT DESCRIPTION 

GENEOT is used in the construction of a depth grid file from bathymetry 
contour files, In this discussion, the term "grid" refers to a matrix of 
depth values whose row and column dimensions represent data point spacing 
in terms of a desired resolution parameter DHN. 
1. Selection of Bathymetry Grid Regions 

GENBOT will create a bathymetry grid file for any region of the geoid 
for which contour data is supplied in the format stated in section A-4, 
The user must select the latitude and longitude of the center of the E 
region of interest, a radius in kilometers and the cell resolution in 


kilometers. GENBOT will attempt to create a square matrix whose semi- 
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axes represent the selected radius and whose elements represent depth 
values in kilometers. 

A 30-km data margin, outside the calculated axis limits, is required 
to begin interpolation. Should this requirement not be met, GENBOT will 
symmetrically reduce the length of one or both axes until the resulting t 


rectangular matrix plus margins with the desired center position intact 


is totally encompassed within the data region. The center coordinates, 

area radius, and cell resolution must be selected so as to avoid generating 

a grid which is impracticably small, In particular, a generated grid may 

not overlap a pole, and one with a pole on a margin boundary must be sup- 

plied data from eighteen contour files to cover 180 degrees of longitude, 
Contour data currently available in usable format cover 40 degrees 

of longitude of the North Atlantic Ocean, specifically the Norwegian and 


Greenland Seas, from longitudes 20W to 20E and latitudes 60N to 90N 


(Fig. 8). The Naval Research Laboratory, Acoustics Division, Arctic Section, 


has additional files which cover the entire North Atlantic land and a 
regions which may be utilized by GENBOT after modification. However, 
GENBOT need not be the exclusive means of grid instruction. The only 
requirement is that the grid exists as a file in the format specified. 
The user may wish to draw on other souces of bathymetry data, such as 
SYNBAPS [8], and tailor other available resources to construct a usable 
grid. 
2. Grid Construction Techniques 

Once supplied with the user's selection for the grid title, center 
coordinates, cell resolution and area radius (TTL, CNLAT, CNLON, DHN, 
ARAD), GENBOT begins by defining a north-oriented square region whose 
boundaries are at a distance ARAD from the center. The subscripts of 
the two dimension matrix depicting this region will represent equal units 
of DHN spacing from the origin along the x and y axes of the square using 
a cartesian coordinate system. Thus, while all distances calculated 
during grid generations are geodetic distances from the center, point 
positions are made relative to the grid origin or southwest corner. This 
technique results in a plane projection which is distortion free along 
straight tracks which pass through the center. Subroutine CORNER then 
calculates and prints a table of the initial cartesian and geographic 
coordinates of the grid corners. The cartesian coordinates will always 
range from (0,0) to their maxium values, initially double the value of 
ARAD. The geographic coordinates should reflect some meridian conver- 
gency. This table will be repeated for the final grid, showing the 
adjustments made, if any, to meet the data requirements. 

Many geographic coordinate pairs in an input file may be well outside 


the area relevant to the grid, so GENBOT calculates gate values called 
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SLAT, ELAT, SLON, and ELON for data discrimination. SLAT and ELAT are i 
always one degree beyond the minimum and maximum latitude values calcu- 
lated by CORNER, but SLON and ELON may extend beyond the minimum and 
maximum longitudes by 1, 10, or 20 degrees, depending upon the maximum 
latitude. 

The grid region and its data environs recognized by GENBOT at this 
stage is shown in Fig. A.1. All points within these environs are sub- 


ject for consideration during interpolation of depth values in the grid. 


FIGURE A.1. Geographic region used for grid construction. 


raves 
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Data files used as input for GENBOT are digitized represenations of 


an NRL bathymetry chart [4] and were recorded on Mobil Oil Company Tape 
NAV-78 for use in a computer program to draw portions of the chart in the 
form of bathymetric contour lines and their depth labels as they exist 

on the chart. The records are in the card image format and consist of 
three types: 

¢ NAME - depth value coordinates 

« CON - contour line header label 

¢ Chart coordinates of points along the contour line 

In the original form the numerical data on the tape consisted of lati- 
tude/longitude pairs, providing no means of determining the depth repre- 
sented by the contour lines without visual inspection of the computer plot. 
The tape used in this research is similar in format to NAV-78 except that 
the contour line header records include depth data (provided by NRL), and 
only the regions shown in Fig. A.2 and listed in Table A.I, are covered. 
The numerals of the data set names were chosen to reflect the original 
file number on the source tape, NAV-78. 

Using the magnetic tape files selected by the user, GENBOT begins 
reading the digitized contour data, stopping after each record to dis-~ 
criminate the points, calculate geodetic distances, transform the data 
to x, y, zcoordinates, and load the coordinates into three axis vectors, 
AX, AY, and AZ (which will be used later for interpolation) while simul- 
taneously keeping a ruining plot of the contour lines. 

After all input files have been processed, the vectors AX and AY are 
sorted into ascending order by the SSP3 subroutine RAIN1 [5]. The mini- 
mum and maximum values are tested to determine if the 30-km outside data 


margin is included in the range. If not, axis reduction begins until 
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25° 20° 15° 10° 5° O° 5° 
FIGURE A.2. North Atlantic regions, with their magnetic tape file numbers 
and data set names, available for use by GENBOT. All are 10° 


by 10°. Other regions available on the same tape are listed 
in Table A.I. 


TABLE A.I 
FILE ONE OF CONTOUR DATA RECORDED ON A MAGNETIC TAPE 


DATA: BATHYMETRIC CONTOURS OF THE NORWEGIAN AND GREENLAND SEAS FROM 
LAT 60N TO 90N and LONG 20W to 20E. 


SOURCE: ADAPTED FROM MOBILE OIL CO. TAPE NO. NAV-78 AND NAVAL RESEARCH 
LABORATORY DATA. 


STRUCTURE: EACH FILE COVERS A 10 X 10 DEGREE REGION. THE FOLLOWING LIST 
IDENTIFIES THESE FILES. THE FIRST VALUE IS THE FILE NUMBER, FOLLOWED BY 
THE DSN, THEN THE LAT/LONG PAIR OF THE SOUTHWEST CORNER (LABOR,LONOR) OF 
THE REGION, AND FINALLY THE NUMBER OF CARD IMAGES CONTAINED IN THE FILE. 


THE TAPE IS STANDARD LABELED, DENSITY 1600 BPI, LRECL=80,BLKSIZE=8C0, 


2 CONF13 +80-20 256 
5 CONF26 +80-20 4343 
8 CONF34 +70+00 12443 
11 CONF35 +70-10 9428 


3 CONF24 +60+10 569 


6 CONF27 +60-10 6001 
9 CONF37 +80+00 491 
12 CONF36 +80+10 409 
4 CONF25 +60+00 3026 
7 CONF30 +70+10 857 
10 CONFO6 +70-20 4454 


13 CONF38 +80-10 972 
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either the margin requirement is met or goes below a 7.5 km semi-axis 
length, the latter case preventing grid construction although a contour 
plot may still be produced. Subroutine CORNER is then called again to 
print the corner coordinate table. If no adjustments were necessary, the 
second table will be identical to the first. 

The interpolation process begins at the southwest corner of the grid 
and proceeds in a columnar direction (i.e., completing the range of x 
values for each Ys value) via successive calls to the SSP3 subroutine 
RAIN2, The accuracy of the interpolated values is proportional to the 
density of data in the given area of interest. Every ten values inter- 
polated are written as a record on the ocutput file. The process continues 
until the total number of matrix elements has been reached, at which time 
the output file is closed and the program is terminated, 
3. User Instructions 

a. JCL and Program Deck. 

The program section of a sample deck is printed below. The TIME 
and REGION.GO parameters were chosen in anticipation of a maximum number 
of 5000 data points in the area of interest, a square region of 210 km 
on a side (150 km plus margin). The program listings included in this 
appendix are part of this example. 

//(Standard job card) , TIME=20 
// EXEC FORTCLGW,REGION.GO=250K (See Appendix 0.) 
// FORT.SYSIN DD * 
(GENBOT, CORNER, GEOPLT, GEODST listings) 
(RAIN1, RAIN3, RAIN2 listings) 
/* 
The following DD statements identify the magnetic tape files 


representing the ten by ten degree regions into which the subroutine 


67 


extends. In the example, the area chosen is centered at coordinates 
70°-04'N/15°-51' (70.07059/15.85302) with a radius of 75 km. This area 
requires two data files as input, CONF24 (file 3) and CONF30 (file 7). 
Since tape data is sequential, the order of the DD cards should be in the 
order of the file numbers (appearing in the LABEL parameter). 
//GO.FTO1LFOO1 DD UNIT=34u0-3, VOL=SER=NPS,DISF=(OLD,PASS), 


//  LABEL=(3,SL,,IN) ,DSN=CONF30 


//GO.FTO1LFOO2 UNIT=AFF=FTOLFOO1, VOL=(,RETAIN, REF=*,GO.FTO1FOO1), 


// DISP=(OLD,PASS) , LABEL=(7,SL, , IN) ,DSN=CONF30 

More information on magnetic tape usage at NPS may be found in Ref. 10. 
Next must follow the DD statements specifying the data set and 

device for permanent storage of the generated grid. The permanent device 
may be another magnetic tape, a disk, or punched cards. In all cases, 
the size of the output data set is important, and the DCB parameters 
should be considered in advance if the grid size is different from the 
example. GENBOT produces card image out put, meaning the logical record 
length is 80 bytes. 

e Punched cards. If the grid is to be stored on punched cards, the 
user need only change the data set reference number on the appro- 
priate WRITE statements in GENBOT from a two to a seven. The num- 
ber of cards punched will equal a 3+[(row dimension) (column 
dimension)/10] rounded up to an integer value. This represents 
three header cards followed by ten depth values per card in columnar 
order, as required by ECTRACE. Thus if the input grid radius was 
75 km and input cell resolution 1 km, the maximum parameter calculated 


by GENBOT will be: 
XMAX 
150 


the 


// 
// 
// 


resulting in a matrix of 22801 elements and a punched output of 

2284 cards. 

Disk. Grid storage on a disk device is the most convenient form for 
multiple ECTRACE runs over the same sea bed. The following DD state- 
ment stores the example grid on a 3330 disk drive device at NPS as a 


data set named $999.GRID. 


GO.FTO2F001 DD UNIT=3330,VOL=SER=DISKO3,DISP=(,KEEP), 

DCB=(FECFM=FB, LRECL=80, BLKSIZE=800) , SPACE=(TRK, 20), 

DSN=S9999.GRID 
Blocked at 800 bytes, the output data requires 14 tracks per 10 
records (card images). More information on the DCB and DSN parameters 
may be found in Ref. 11. When the generated grid is no longer needed, 
and if the user wishes to generate a new grid, or if space has already 
been allocated on the device, the DISP parameter should be SHR in 
place of (KEEP) and the DCB and SPACE specifications may be omitted. 
Disk data sets are periodically purged by center operations persbdnnel, 
but obsolete data set space should be cleared or reallocated by the 
user, 

b. Input Data, 


The following input data are user selected values and follow 


//GO.SYSIN DD* card. 
Card Variable FORTRAN Field Descriptors Remarks 
1-2 TTL 6A8/6A8 Title 
3 CNLAT, 4F10.5 center latitude, 
CNLON, center longitude, 
DHN, ARAD cell resolution, 
radius 
4 NFILES, LEN 215 : Number of input 
files, plot 
length 
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The magnitude of LEN determines the plot size (<19 inches). If 
LEN is negative, no grid will be constructed. If its magnitude is less 
than four, no plot will be produced. 
ec. Printed Output 
GENBOT printed output supplies the following: 

- A repeat of the user input data, 

« A table of the grid corner coordinates and the depth matrix 
dimensions using the input parameters. 

- A summary of results for each contour line, including its depth, 
and the title from the header label, the number of points on the 
line (NPTS) and the number of points used from the line (KPTS). 

- A corner coordinate table of the final grid. 

- A summary of the interpolation results. A line is printed for 
each interpolated point whose accuracy may be questioned. The 
meaning of the IER parameter may be found in the RAIN] listing 
in Appendix I. 

d. Plotted Output 
The contour lines of the input file are plotted with the local 
parallels and meridians overlaid. The resultant subregion represented by 


the output grid is outlined by a rectangle with labeled axes. The axis 
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values are referenced to the desired grid origin and not the resultant 
grid origin, the latter being at (0,0) always. 
e. Programming Notes on Contour Data Interpolation 
Subroutine RAIN2, the interpolation routine used by GENBOT, is 
aptly suited to interpolate points that are more or less randomly spread 
in both horizontal directions. The nature of contour data, however, is 


1’ 


typically a collection of closely spaced depth values along widely spaced 
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contour (compare Figures III.1.-and-A.l.). This type of data arrangement 
frequently results in a colinearity problem for RAIN2, which can be resolved 
only by searching ever widening data cells for a non-colinear pcint, result- 
ing in a degradation of accuracy. In addition, the routine is highly sensi- 
tive to contour line kinks which result in right-angle deviations in the 
interpolated contour lines (which may be plotted by GRDCHK). 

Correcting the deficiencies of an interpolating routine is beyond 
the scope of this research, but the Naval Research Laboratory, Arctic Sec- 
tion, is currently investigating methods for generating grid-matrices from 
contour data. Users are encouraged to modify the code used in GENBOT as 
desired to produce results which more closely reflect the actual data 
supplied. 

4, Input Contour Data Format 

As stated previously, some contour data already exists in a format 
usable by GENBOT, but if data is supplied by another source, it must be 
arranged in a file in the manner and format specified below. All records 
are 80 byte card images. 

The first record of the file contains data on the file itself. GENBOT 
uses a 415 field to represent the origin latitude/longitude pair, the num- 
ber of contour lines in the file, and the number of card images it may 
expect to read. 

Contour line data exists as separate sequences of cards and may 
appear in any part of the file as long as each line or line segment is 
in an unbroken sequence with a header card, The header card must contain 


the line's depth in hundreds of meters, right justified in a 14 field in 


columns 41-44. The cards following the header card contain the latitude/ 


longitude pairs in degrees and decimals of degrees of the points along 
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the contour line in the form of four pairs of coordinates (in FORTRAN 
fields 8F10.5) per card. South and west values are negative. A value 
greater than 90.0 in magnitude in any latitude field is interpreted as the 
end of data for the contour line, and GENBOT expects the next record to 

be either another contour line header label, the end of the file, or a 
depth value position record, 

Depth value position records are optional and are used during the 
chart-drawing phase of the program to place alphanumeric symbols on the 
chart product. They may exist any place after the first record in the 
file except imbedded in a contour line sequence. Columns 4-8 must con- 
tain "NAME" in an A4 field, columns 19-38 contain the latitude and longi- 
tude of the symbols in a 2F10.5 field and the symbols themselves must be 
in columns 46-47. These records are normally used to draw depth values 
in hundreds of meters on the contour lines, 

If the user wishes to use GENBOT on externally supplied contour data, 
it it-recommended that a dump be performed on one of the smaller files 


(such as CONF13) of the GENBOT data tape to verify the required format 


for adaptation. 
5. Qutput Grid Data Format 

GENBOT constructs the output file in accordance with GRDSCT require- 
ments, Since GRDSCT is the subroutine used by TRACER (of ECTRACE), 
G3DP and GRDCHK, any externally supplied grid data must all be of the 
same structure as outlined below in order to be used by those programs. 
Again, all records are card images. 

a. The first two records contain the input grid title in the 10A8 


fields. 


: b. The third record contains the input latitude and lontitude of the 


grid center, the point spacing, and the number of rows and columns of the 
matrix representing the grid. The record format is (3F10.5,215). 

c. The fourth and remaining records contain ten depth values per 
card, in kilometers, as 10F8.4 fields, beginning with column one (y = 0) 
and proceeding all x-values before going to the next column, in accordance 


with the FORTRAN convention of columnar storage of matrices. 
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APPENDIX B 
SYNGEN DESCRIPTION AND USER INSTRUCTIONS 


1. Program Description 


a. SYNGEN is a collection of simple functions which may be selectively 
used to generate and record a matrix of depth values representing a desired 
synthetic bottom configuration. Six bottom cases are provided (see below), 
but the program may be easily modified by the user to produce a case not 
described here. A program listing is provided in Appendix H. 

The user selected parameters are: 
¢ CNLAT - any desired center latitude value. 
CNLON - any desired center longitude value. These parameters 
are requested in compliance with the required output file 


format. 


DMIN ~ the minimum depth value (D,) to be attained within a 


75 km radius from the grid center. 
DMAX - the maximum depth value (D,) within the same radius. 
ITYPE - an integer code for the bottom type desired. 

b. Table B.1 lists the ITYPE code, the bottom type, the functions 
used to generate the depth matrix, and the matrix dimensions. All grids 
represent a region covering a 150 km in both the x and y directions, but 
the matrix size will be 11 by 11 or 151 by 151, depending upon the bottom 


type (DHN values are 15 km and 1 km, respectively). In all cases, 


c. The depth file provided by SYNGEN conforms to the format required 
by ECTRACE, GRDCHK, and G3DP. When one of these programs uses a SYNGEN 
product, the array dimension in the program's source code must be checked 


to conform to the dimensions of the generated matrix. 
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Conical 
Seamount 


Sinusoidal 
Undulation 


Parabolic 
basin 


TABLE B.I 
SYNTHETIC SEA BED FUNCTIONS 


Matrix 


Function Parameters Dimensions 


- ¥ y = ADd/135 11 x 11 
1 


Os xs 135 


+] 6 (x) | y = AD/15 
elsewhere 6(x) = x - 75 


60 s< x < 90 


-|y6(x) | 6 = AD/15 
elsewhwere 6(x) = x - 75 


60 < x § 90 


D + yr y = AD/75 


Do elsewhere r = kx ~ 75)%+ (y =75)* 


D ylitcos{5(x) /R}] y = Ap/2 


6(x) = 2(x - 75) 


R= 15 


y = AD/(75)? 151 


r= ‘e - 75)*+ (y - 75)? 


2. User Instructions 
THE JCL of a sample deck are arranged as follows: 
//(Standard JOB card) ,TIME=2 
// EXEC FORTCLG, REGION.GO=250K 
//FORT.SYSIN DD * 
(SYNGEN source code) 
/* 
//GO.FTO2FOO1 DD UNIT=3330,VOL=SER=DISKO3, 


// DSN=S9999.SYNBED, DCB=(RECFM=FB ,LRECL=80 , BLKSIZE=800) , 


| 
. 
// DSP=(NEW,KEEP), SPACE=(TRK, 20) 
//GO.SYSIN DD * 
(user selected data) 
is } 
Notice that the output data set is directed to DISKO3. The DSN is the | 
user's selection for the data set name in accordance with NPS system re- | 
quirements. A 151 x 151 depth matrix requires twenty tracks of disk i 
space on a 3330 device. | 
Only one data card is needed to supply the program with the user's 


parameter selections. Its field descriptor is (4F10.5,15) and contains 


the values reflected in the same order as listed in paragraph l.a. 
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APPENDIX C 
G3DP DESCRIPTION AND USER INSTRUCTIONS 


1. Program Description 
a. G3DP uses the NPS library subroutine PLT3Dl1 [Ref. 12] to draw a 


perspective surface plot of a subset of a depth matrix. This plot is a 
normal projection of the surface using a focal length and line-of-sight 
calculated from the user's selections for the viewpoint position angles 
DEGUP and DEGCW. As the variable names imply, the viewer is elevated 
DEGUP degrees from a reference depth and rotated DEGCW clockwise from a 
south-to-north viewing aspect. For example, the standard working matrix 
(dimensioned 151 by 151) is best viewed with DEGUP and DEGCW values of 
30 degrees and 45 degrees respectively, providing the same viewing aspect 
demonstrated in all the the G3DP products printed in this report. Using 
a DEGUP value of zero results in an isometric rather than a perspective 
projection. The DEGCW selection should be between zero and 89 degrees 
(south-north to west-east viewing aspects). 

b. As in ECTRACE and GRDCHK, G3DP requires source code modification 
by the user to tailor the program to the desired dimensions of the work~ 
ing matrix. Instructions in this regard are included as comments in the 
FORTRAN listing (Appendix K). 

ec. An important characteristic of G3DP is that it uses only the 
odd-numbered rows and columns of the working matrix to produce the sur- 
face plot. This reduction in resolution is required to reduce core re- 
quirements and unclutter the graphic output caused by typically large 


(greater than 75 by 75) working matrices, 
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2. User Instructions 


The JCL of a sample program deck are arranged as follows: 


//(Standard JOB card) ,TIME=2 


//EXEC FORTCLGW,REGION.GO=300K (Appendix 0.) 


//FORT.SYSIN DD* 


/* 


(G3DP source code) 


//LINK.USDD UNIT=3330, VOL=SER=DISKO1,DSN=1570. ECTRACE, 


// 


DISP=SHR,LABEL=(, , ,IN) 


//LINK.SYSIN DD* 


/* 


INCLUDE USDD(GRDSCT) 


ENTRY MAIN 


//GO,FTO1F 


//GO.SYSIN DD* 


/* 


(user-selected data) 


The LINK step provides access to GRDSCT for loading the working 


matrix 


from the input file. The input data set is the same example used 


in ECTRACE and GRDCHK. The card input data are user selected and con- 


tained 


on one card, as follows (FORTRAN field descriptor 5F8.3) 

XST - working grid matrix x-origin relative to that of the 

input grid (km). 

YST - working grid-matrix y-origin relative to that of the input 
grid (km). 


DEGUP - viewpoint elevation angle (0 § DEGUP < 89 degrees). 


* DEGCW - viewpoint rotation angle from south to north 


(O < DEGCW < 89 degrees). 


¢ LW - width of the plot (5 < LW < 18 inches). 


APPENDIX D 
GRDCHK-USER INSTRUCTIONS AND OUTPUT DESCRIPTION 


pee Purpose and Output Description 


Cases may arise where the accuracy of depth values in a generated ma- 
trix may be in question because of data sparsity, interpolation deficiencies 
or other reasons. Program GRDCHK assists in evaluating the generated grid 
by printing and plotting selected columns of the matrix for comparison with 
that of the source contours. The contour plot may also be used on a light 
table with an ECTRACE or ECCOM horizontal-ray trace plot (Figs. 5, 6, and 
7). 


The NPS library routine CONTUR is called for producing the contour 


plot. Subroutine GEOPLR is a modified version of GEOPLT to accommodate 
the rotated plot produced by CONTUR.. | 
2. User Instructions 
To aid in this discussion refer to a computer listing of GRDCHK code 
in Appendix J. A listing of the CONTUR subroutine may be found in Ref. 6. 
The following example applies to a standard 151 by 151 working grid. 
Below is a listing of the JCL cards and the position of the program deck. | 
//(Standard job card), TIME=5 | 
// EXEC FORTCGLW, REGION.GO=350K | 
//FORT.SYSIN DD * | 
(GRDCHK ,CORNER,GEOPLT,GEODST source decks) 
/* 
The following DD statements link to a member of ECTRACE for allocating 
a working subregion of the input. 


//LINK.USDD DD UNIT=3330,VOL=SER=DISKO1, 
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// DSN= 1570, ECTRACE, DISP=SHR,LABEL=(,,, IN) 
//LINK,.SYSIN DD * 
INCLUDE USDD (GRDSCT) 
ENTRY MAIN 
/* 
Next follows the DD statement identifying the depth matrix input file. 
The same example is used for the GENBOT output file. 
//GO.FTOFOO1 DD DSN=S9999.GRID, UNIT=3330, VOL=SER=DISKO3, 
//  DISP=SHR, LABEL=(,,,IN) 
Finally the user input parameters are included. The card 
//GO.SYSIN DD * 


is followed by: 


Card Parameter Names FORTRAN field description 
i NL, LEN 215 
2 ZST, ZEND, XMIN, 4F6.1 
YMIN 
3 IPS, IPE, JPS, JPE 415 


where 

* NL is the number of contour levels desired. The shallowest and 
deepest levels, if known, should be read into ZST and ZEND. If 
the boundary levels are not known, NL should be made negative, 
and subroutine CONTUR will calculate the boundary and contour 
levels. If NL is greater in magnitude than 100, no contour plot 
will be produced. 

¢ LEN is the length of the contour plot's longest side in inches. 
If LEN is negative, a cartesian grid will he drawn on the con- 
tour plot instead of geographic coordinate lines. If LEN is 


less than 4, no contour plot will be produced. The maximum 
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magnitude of LEN is 19 inches. 
ZST AND ZEND are estimates of the smallest and largest values in 
the working matrix (ZB). 


XMIN and YMIN establish the origin of the working matrix relative 


oat a Lim teahinabai ae: 


to that of the input matrix. The program assignments for KMX and 
KMY, the recorded cell spacing DHN, and the dimensions of the work- 
ing matrix determine the subregion size parameters in the same man- 
ner as they do in ECTRACE, since subroutine GRDSCT is linked to 
this program. In effect, the user must tailor certain portions 

of the program cards to compensate for the FORTRAN G inability 

to perform object time dimensioning. 

IPS, IPE, JPS and JPE identify the beginning and end rows and 
columns of the working depth matrix to be singled out for inspec- 
tion. Thus, if ZB is dimensioned (150,150) and the above values 
are 1, 150, 74, and 78 respectively, GRDCHK will print the depth 
values of these 750 elements of the ZB matrix and plot five 
columnar profiles (Fig. 4). If this record is omitted, only 


the contour plot will be produced. 
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APPENDIX E 
ECCOM OUTPUT DESCRIPTION AND USER INSTRUCTIONS 


l. Description of Output 


ECCOM allows the user to combine several selected ECTRACE runs and 
eliminate specific rays as desired to produce an improved horizontal ray 
path plot for investigation of curvature (Fig. 7 and Fig. 18). The size 
of the ECCOM plot will be the same size as first ECTRACE run of the input 
card deck unless the user opts to alter the scale factor (FACT) to control 
the plot dimensions. 

The printed output is a list of all the rays processed for the ECCOM 
plot, giving each ray's initial elevation angle and heading. Additionally, 
the number of points used to plot each ray path is printed. This number 
of points corresponds to the number of reflections plus the starting point 
of the ray. Refer to Table VII for a sample listing of data for the ECCOM 
51 plot, shown in Fig. 7. 

2. User Instructions 
The following JCL cards and parameters are required for an ECCOM run: 
¢ JOB card ~ The CPU time parameter should allow four seconds for 
every 55 rays to be plotted. 
EXEC card - This card must specify FORTCLGW procedure. 
Calling program - Following the 
//FORT.SYSIN DD * 
is the ECCOM program deck. See Appendix L for the program listing. 
Data cards--Following the 
//GO.SYSIN DD * 
is the sequence of punched output decks, each individually arranged 


as described in Chapter V, Section C. The user may change the title 
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card (card #1) of the first ECTRACE punched deck to identify the 
ECCOM run title. In addition, the user may change the scale factor 
(FACT) on card number 4 to modify the plot size. FACT was calculated 
by ECTRACE in response to the user's selection for the horizontal 
plot length (Chapter V, Section B). The following are sample values 


of FACT and recommended limits: 


FACT Horizontal Plot Axis Length (inches) 

0.333 5 (minimum recommended) 
1.000 15 

1.267 19 (maximum recommended) 


If the user desires a plot size outside the range recommended, re~ 


fer to Ref. 7. 


Cards 1 through 4 of each ECTRACE deck after the first are ignored 
by ECCOM, but the user should verify that they all reference the same 
portion of the same input grid-matrix in order for the ECCOM plot to be 


meaningful. 
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APPENDIX F 
ANALYTICAL VERIFICATION OF ECTRACE HORIZONTAL RAY PATHS 

The accumulation of ray heading changes after repeated specular 
reflection from a sloping sea bed describes a faceted horizontal path 
whose apparent curvature may in some instances be analytically described. 
Harrision [2] derived some solutions for horizontal ray paths over simple 
bottom topography cases. Most of these cases may be modeled by SYNGEN and 
used to test the validity of the horizontal ray paths mapped by ECTRACE. 
We chose an isovelocity sound speed condition and a conical seamount because 
its surface geometry effectively challenges the facet methed of topograph- 
ical modeling. 

For a conical seamount of arbitrary slope but whose apex just touches 
the sea surface, Harrison gives the ray paths in polar coordinates (r, %) 


as 


45 


r=r (1 - cos*@ cos*¢ ) 
0 0 0 


(1 - cos?@,cos%4) @ 


* cosec (o + o') 


sing,cos8, 


initial distance from the apex 
initial deviation angle 
initial azimuthal angle 


a constant given by setting $= 0 when r = rT 


sind,cos6, sin"}(1 - cos*@ cos*> Ne 


(1 - cos’9 cos’g,) (F.2) 
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Figure F.l.is a diagram of a representative ray path. a is the 


deflection angle and represents the asymptotic angle of the shadow zone of 


the ray fan, 


| 
Mi 
ae 
Kt 
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FIGURE F.1. Polar representation of the horizontal projection of a ray 
path for a conical seamount. 
In the test run (Trial 83, Fig. F.2), the source was positioned at 
a range r, of 50 km from the apex, bearing 180°, at a depth just below 
the surface. The receiver was placed at the apex so that for each ray 
the ECTRACE range TE could be read directly from the printed ray history 
table values for DIST TO REC and ¢ is the reciprocal of BNG TO REC. A 


zero gradient was used. 
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Two ray lines were projected using 8.5 5 degrees for ray fan 1 and 
@, = 10 degrees for ray fan 2. > ranged from 270° for the first ray 
to 90° for the last ray in each fan, with an increment of 15° between rays. 
Two rays were chosen for analysis, ray 6 (6, = 5 degrees, $,= -15 degrees) 
and ray 21 (6, = 10 degrees, $, = 15 degrees). The relevant portions of 
their trace histories are reproduced in Table F.I and Table F.II. Only 


the surface reflection values are used, where the slant range DIST becomes 


the horizontal range r 


Inserting 6 and @ values into Eq. (F.1), 
° ° 


E° 
loaded in a programmable calculator, we obtained a table of values for r 


for both rays. The results are shown in Table F.III and Table F.IV along- 


side the ECTRACE values TE and the observed error € . These same 


results are plotted in Figures F.3 and F.4 (rvs 9%). 


TABLE F.I 


RAY 6 (6 OF RAY TAN 1) TRACE HISTORY 


REV DEPTH GRID ELEV BRNG DIST 
TYPES HEADING ANGLE TO REC FM REC 


344.500 
344.500 
343.010 
343.010 
340.647 
340.847 
338.231 
338.231 
336.557 
336.557 
335.915 


335.915 


TABLE F,II 


RAY 21 (8 OF RAY FAN 2) TRACE HISTORY 


REV. DEPTH GRID ELEV BRNG DIST 
NO. TYPES (KM) HEADING ANGLE TO REC FM REC 
0 INIT 0.001 15.000 10.000 360.000 50.000 
1 BOTT 2.451 15.649 -17,119 354.386 36.838 
2 SURF 0.0 15.649 17.119 348.770 29.483 
3 BOTT 1.654 17.253 -23.468 343.153 24.866 
4 SURF 0.0 17.253 23.468 337.518 21.761 
5 BOTT 1.303 20.043 ~28,702 331.883 19.590 
6 SURF 0.0 20.043 28.702 326.243 18.046 
7 BOTT 1.128 24.047 -32.594 320.603 16.949 
8 SURF 0.0 24.047 32.594 315.014 16,200 
9 BOTT 1,046 28.867 ~33.973 309.425 15.723 
10 SURF 0.0 28.867 33.973 303.767 15.480 
11 BOTT 1.029 34,081 ~33.685 298.105 15.459 
12 SURF 0.0 34,081 33.685 293.479 15.661 
13 BOTT 1.071 38.566 -31.042 286.852 16.097 
14 SURF 0.0 38.566 31.042 281.203 16.802 
15 BOTT 1.187 41.986 ~26,620 275.55 17.834 
16 SURF 0.0 41.986 26,620 269.889 19,295 
17 BOTT 1.420 44,169 ~20.799 264.224 21.345 
18 SURF 0.0 44.169 29.799 258.538 24,279 
19 BOTT 1,904 45.278 ~14 066 252.851 28.620 


20 SURF 0.0 45.278 14.066 247.160 35.467 


$ 
0,0° 


-22.736 
-45.408 
-68.425 
-91.462 
-114,49 


-137.46 


TABLE F.III 
Trial 83 Ray 6 


(® = $, constant) 


TABLE F,.IV 


Trial 83 Ray 21 


$ r Tr € 

0.0° 50.0 km. 50.0 km 0.0 km 

14.230 29.470 29.483 -0.013 

22.482 21.745 21.761 -0.016 

33.757 18.026 18.046 -0.020 | 
44.986 16.178 16.200 -0.022 | 
56.233 15.548 15.480 -0.022 | 
67.521 15.643 15,651 -0.018 
78.797 16.791 16.802 -0.011 
gO.111 19.309 19.295 0.014 

101.46 24.365 24.279 0.086 

112.84 35.822 35.467 0.355 

124.219 75.857 73.529 2.238 

Ray 20 ™ * 8.687 km (% = 0, constant) 
ce 
v4 * "Eee = -0.005 km 
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Harrison verificaton plot of range vs % for ray 6, 


FIGURE F.3. 
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FIGURE F.4. 
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Harrison Verification Plot of range vs $ for ray 21 
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These results show very close agreement up to a distance well past 
the closest point of approach, and only begin to diverge slightly when the 
rays approach their asymptotic values for %. Investigation of other rays 
of the same trial supports this observation. Clearly the fact that the 
plane facets do not all have the same slope as the cone they collectively 
represent, but vary as much as 4% between their extremes, explains most 
of the error. The validity of the model decreases near the apex, where 
fewer data points and hence plane facets are used to approximate the 
conical surface. In addition, the progressively shallower grazing angles 
assumed by the rays after CPA cause the horizontal paths to become more 
sensitive to this facet effect. 

Reference 13 also provides an equation for the closest approach of 
the ray fan envelopeto the seamount center, developed by letting 


$, = 0 in Eq. (F.1). Thus 


y 
it 


r. sin 6 
ce q 4 


8.682 km. (F.3) 


Ray 21 of Trial 83 was the radial ray (>, = 0°) and its CPA on 
the printout was 8.687 km, a difference of only four meters after il 
reflections from the bottom model. The ray should have reversed its 
heading exactly 180 degrees after CPA, but actually experienced a 174 
degree heading change, revealing the sensitivity of the model to radially- 
directed ray paths. 

Lastly, Ref. 13 defives the asymptotic angle of the shadow zone of 
the ray fan (Fig. F.1). A plot of the final azimuthal angle >, versus 
8, for ray fan 2 (Fig. F.5) of the test run exhibits an interpolated 


value for a of 40 degrees for >, = 29 degrees. Testing this value with 
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FIGURE F.5. 


the formula 


1 
a # 3(nsin’e,) /s 


= 39,187 degrees (F.4) 


shows close agreement with the theoretical value. Figure F.6, an ECCOM 
composite for the seamount, is a graphic example of the dependence of 
@ on the initial elevation angle. 

It is worth noting that the analytical expressions for horizontal 
projections of ray paths apply only to exact functional bottom topographies 
and that the curved three dimensional ray path surfaces only approximate 


an actual surface of vertical planes joined at bottom reflection points. 
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FIGURE F.6. ECCOM Trial 28 showing the effect of a synthetic conical 
seamount on ray propagation. The ray's initial headings 
varied from 20° to 58° and the initial elevation varied 
from 4° to 18°. Selected rays from five ECTRACE runs 
(100 rays) were combined. 
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APPENDIX G 
TABLES 
TABLE I 
A SUMMARY OF THE PROGRAMS AND THEIR COMPONENT SUBPROGRAMS 


GENBOT - plots source contour data and generates depth matrix 
CORNER - calculates boundary coordinates 
GEOPLT - plots geometric grid lines 
GEODST - solves geodetic triangle 
RAIN 1/2/3 - perform sorting and interpolation (from SSP3) 


SYNGEN - produces synthetic sea bed depth matrices 


G3DP - produces perspective surface plot of depth matrix 
A. GRDSCT ~ loads working surface plot of depth matrix 
B, PLT3D1 - draws surface plot (from NPS) 


GRDCHK - allows inspection of depth matrix and draws contour plot 
of same 
GRDSCT = see above 
CORNER - see above 
GEOPLR -— same as GEOPLT but for a rotated plot 
GEODST - see above 
CONTUR - draws a contour plot of depth matrix (from NPS) 


ECTRACE - a collection of subprograms which perform ray tracing 

A. TRACER - controls all subprograms and manages ray traces 

GRDSCT - see above 

IDTSUB - identifies subscripts of depth matrix 

DEEP - calculates depth at given coordinates 

IDPROF - identifies water mass 

NUPROF - calculates water mass parameters 

CHNLIM - calculates sound channel limits 

CONTAC - calculates point of bottom contact 

BOANG ~ calculates direction and grazing angle of bottom 
reflected ray 

FBLOSS - calculates bottom loss 

ANGPRT - prints ray fan bottom bounce table 

BNGPLT - initializes plotting and draws plot borders 

BDTPLT - draws bathymetry profile 

RNGPLT - draws vertical plot 

T2DPLT - draws horizontal plot 

SSPPLT - draws sound speed profiles 

ENDPLT - draws plot legend and terminates plotting 
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TABLE II 


INPUT PARAMETERS TABLEAU (COMPUTER PRINT OUTPUT) 


SOUND SPEED PROFILES 


vAQ) 


0.0 
50.00 
60.00 
70.00 
90.00 

500.00 
1050.00 
1500.00 
4000.00 


ZC 


0.0 

50.00 
200.00 
375.00 
750.00 
975.00 
1250.00 
1500.00 
4000.00 


zc 


0.0 
1500.00 
4000.00 


NWMR (Northwest Mohns 


SEMR (Southeast Mohns 


ECTRACE TRIAL 67 *** SYNTHETIC BASIN *** 


C(ZCc) 


1461.08 
1461.00 
1447.00 
1445.00 
1447.00 
1454.30 
1461.00 
1470.00 
1511.70 


C(ZzCc) 


1490.00 
1482.00 
1474.50 
1474.00 
1466.00 
1465.00 
1467.00 
1470.00 
1511.70 


C(ZC) 


1470.00 
1470.00 
1511.70 


ISsOvV (Isovelocity) 


(Table II continued) 
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Ridge) 


Ridge) 


GRAD 


0.0 
-1.40000 
-0.20000 

0.10000 

0.01780 

0.01218 

0.02000 

0.01668 

0.01668 


GRAD 


-0.16000 
~0.05000 
-~0.00286 
0.02133 
-0.00444 
0.00727 
0.01200 
0.01668 
0.01668 


GRAD 


0.0 
0.01668 
0.01668 


Pi pee oe 
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continued TABLE IT 


FRONT/BOUNDARY INTERSECTIONS 
FRONT FRONT 
ENTRY EXIT ENTRY EXIT 
( 0.0, 50.0) (100.0,150.0) ( 50.0, 0.0) (150.0, 100.0) 


GRID FILE CALLED £&&44AAAKARAAKAKRARARARRARRRR RARER RII I IRE RI IR IRAE 

SYNTHETIC SEA BED 

COMPUTER GENERATED BASIN 
RAKKRAKARERERER ARERR EREAERRRRE AREA ARRAERERERRRERREER 

HAD BEEN READ. CENTER LAT/LONG IS 70.07059/ 15.85032 AT 

(75.000, 75.000) KM. RANGES ARE 0.0 TO 150.0 KM IN Y (22801 ITEMS). 


he 


* 
* 
x 
* 
* 
* 
* 


THE FOLLOWING ADDITIONAL PARAMETERS WERE INPUT: 


SOURCE DATA: FREQUENCY 150.00 HZ,GRID COORDINATES ( 40.00,40.00), 

DEPTH 0.850 KM. 8 RAY FANS WILL BE TRACED, FROM -30.0 TO 30.0 DEGREES 
ELEVATION, INCREMENTING EVERY 8.6 DEG. 1 RAYS MAKE A RAY FAN FROM 10.0 
TO 10.0 DEGREES TRUE HEADING, INCREMENTING EVERY 0.0 DEGREES. 


RECEIVER DATA: GRID COORDINATES ( 125.0, 40.0), DEPTH 1.749 KM. 


GRID ARRAY: MAXIMUM ARRAY SIZE (151,151).DEPTH VALUES RESOLVED BY 1.0 
KM. SPACING. 


TRACE PARAMETERS: PROPAGATION WILL CONTINUE UNTIL A BOTTOM BOUNCE LOSS 
OF 150.0 DB OR 21 BOTTOM CONTACTS EXCEEDED FOR EACH 
RAY. PATH INCREMENT IS 100.0 M. POSITION ERROR AFTER 
BOTTOM CONTACT WITHIN 0.05000 KM. CHANNEL RANGE 
LIMIT 200.0 KM. SMALL ANGLE 2.0 DEG. 


PLOT OPTION SELECTED. 


PLOT PARAMETERS: x Y vA R 
AXES: 0.0 0.0 -1.0 -10,0 
10.0 10.0 2.0 7.00 
150.0 150.0 6.0 123.00 


SCALING: LENGTH 4 IN., FACTOR=0. 2667 


EACH SEGMENT REPRESENTS 10 POSITION CALCULATIONS. 


(program parameters are user supplied via DATA cards) 
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TABLE III 


Ray History Tableau ECTRACE Trail 67 (computer print output) 

Ray 8 (of Ray Fan 8) TRACE HISTORY 

REV. NO. TYPE Y (KM) DEPTH (KM) HEADING ELEV. ANG. 
10.000 30.000 
12.533 -20.323 
12.533 21.104 

04 14.157 -14.658 

0 14.157 17.615 

34 15.556 -14.938 

0 15.556 18.425 
7 
0 


°o 


62 17.472 -22.313 
17.472 26.133 


8 20.496 -31.868 
20.496 33.698 


6 
116.479 0 
122.046 

126.024 

129.119 

131.279 

133.064 

134.260 

135.305 

135.934 

136.511 

136.749 

136.977 

136.890 

136.802 

136.378 

135.929 

135.080 

134.129 

132.631 

130.845 

128.228 

125.020 

120.443 

114.414 

106.392 

96.179 

84.615 

103.807 73.002 
107.968 62.549 
111.257 954.284 
113.497 47.940 


19 25.415 41.798 
0 25.415 42.680 
82 33.460 ~51.533 
0 33.460 51.978 
6 46.953 -60.129 

-0 46.953 60.372 
4 
0 
4 
0 
5 
7 


e 


69.313 -65.773 
69.313 65.930 


97.778 -66.103 


° 


97.778 66.244 


e 


121.284 -61.160 
121.284 61.347 


e 


135.723 -53.031 


: 144.372 -43.592 
144.372 44.152 


149.526 -34.065 


e 


0 
0 
el 


149.526 35.164 


. 153.024 -26.290 
. 153.024 26.244 
155.238 -19.001 


155.238 20.346 


e 


156.763 -15.037 


° 


4 
0 135.723 53.338 
0 
0 


0 
6 
0 
4 
2 
. 156.763 18.321 
2 158.297 -17.529 
0 158.297 20.418 
48 160.557 25.605 
0 160.557 28.126 
115.211 43.084 713 164.103 35.243 
116.231 39.501 -0 164.103 36.331 
117.049 36.630 2.188 169.812 -45.355 


BOTTOM LOSS GATE EXCEED. TERMINATING RAYTRACE AFTER 3311 
SEGMENTS. PLOT BASED ON 367 POINTS, 41 REVERSALS AND 21 

BOTTOM BOUNCES. MINIMUM DISTANCE TO RECEIVER, 8.57 KM, 

OCCURS AT TIME 166.244 


OSONOWOFOFOWONONOKFOF OF OF OFORP ORF ONONOWOLOLO 


(TABLE III continued) 
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continued TABLE III 
BRING DIST WTR DEPTH... LIMITS 
MASS ______LYR 


en RCE ae ROE 2 Le ee 
0.0 0.0 0.0 
190.00 4.759 3.270 
191.646 11.906 8.802 
192.056 23.308 16.265 
192.809 34.416 24.931 
193.189 48.043 34.520 
193.690 60.783 43.604 
193.961 71.155 51.036 
194.310 78.932 56.802 
194.572 84.778 61.279 
194,812 88.984 64.648 
195.015 92.279 67.440 
195.276 94.626 69.600 
195.483 96.573 71.568 
195.741 97.936 73.103 
195.962 99.130 74.624 
196.234 99.921 75.814 
196.480 100.647 77.058 
196.784 101.058 78.109 
197.074 101.451 79.258 
197.435 101.559 80.307 
197.796 101.669 81.494 
198.245 101,483 82.550 
198.721 101.290 83.914 
199.322 100.758 85.183 
200.003 100.176 86.827 
200.882 99.146 88.467 
201.954 97.956 90.673 
203. 388 96.130 93.144 
205.222 93.996 96.346 
207.779 90.926 100.139 
211.412 87.232 105.373 
216.510 82.606 111.713 
223.790 77.898 110.854 
232.820 73.832 128.722 
242.651 71.917) 137.727 
251.646 71.615 145,718 
258.664 72.722 152.236 
263.834 73.929 157.358 
267.652 75.297 =161.318 
270.375 76.237 164.413 
272.504 77.134 166.955 


kaKKK 
kaKKK 
kkRKK 
7.093 
7.093 
7.526 


7.526 
KEKKK 


keKKK 
kkk 
keke 
kkk 
kkkkKK 
kkekkk 
kkekekk 
kkekkk 
RAKE 
kaekKK 
kkk 
kk 
kkKKK 
kkkkk 
keRKK 
RAKKE 
ReRKK 
kRKRKK 
RkKKKK 
keEKK 
kkkkk 
keEKKK 
REKKK 
kkREKK 
kkkkk 


7.469 
7.469 
8.687 


8.687 
kaRKK 


REKKK 
RREKK 
RRRKK 
RRERK 


e eee #8 @ oe 


qoooooo0oo0coococcoccooceacaococeocoooecoe°dcoco0coococoo0oo0ooceo0co0o00o00°o 


ee © © e «© @# @ @ 


NPN HOrF OF OF Or Or OFrORr ORY RNR SR Or OF Or Or Or Or OF OU 
ogoooo0oocoococeoocceoocooedcooc°ocoe°coocococococo0ocooooo°%o 


*For this run range and bearing of the rays were referenced to the 
source. As explained in Chapter V, in the current version of 
ECTRACE, these values are now referenced from the receiver (REC). 


(TABLE III continued) 
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continued TABLE III 


c GRAZ BIM BTM 

(KM/S) ANGLE ANGLE LOSS LPS 
1.46556 
1.49770 24.081 4.678 4.1 0 
1.49000 
1.51242 16.731 3.392 6.1 0 
1.49000 
1.51743 14.465 2.756 7.2 0 
1.49000 
1.50549 19.479 3.935 10.0 1 
1.46100 
1.49135 27.785 5.015 15.1 1 
1.46100 
1.48151 37.224 5.651 22.2 1 
1.46100 
1.47539 46.949 6.015 30.0 0 
1.46100 
1.47177 56.159 6.182 40.9 0 
1.46100 
1.46980 63.522 6.272 50.9 0 
1.46100 « 
1.46904 66.735 6.262 60.0 0 
1.46100 
1.46963 64.227 6.248 70.9 0 
1.46100 
1.47145 57.381 6.141 80.9 0 
1.46100 
1.47477 48.418 6.040 89.9 0 
1.46100 
1.48045 38.786 5.639 97.3 0 
1.46100 
1.48941 30.571 5.162 103.0 2 
1.49000 
1.50255 22.176 4.155 106.6 0 
1.49000 
1.51582 16.287 2.836 108.3 0 
1.49000 
1.51606 16.298 2.906 110.1 1 
1.49000 
1.50313 22.375 4.256 113.8 0 
1.47000 
1.49023 30.985 5.212 119.5 0 
1.47000 


1.48147 40.627 5.821 127.3 0 


TABLE IV 


Ray fan Summary Tableau for ECTRACE Trail 67 ray number 8. (computer print 
output) 


ANGLES BEFORE BOUNCE. INITIAL ELEVATION ANGLE = 30.00 


HEADINGS : 10.00 
BOUNCE 
27.78 
18.77 
13.96 
16.60 
23.03 
32.51 
42.10 
51.66 
60.18 
65.80 
66.12 
61.15 
52,96 
43.36 
34.76 
25.30 
17.50 
15.03 
19.09 
26.64 
35.75 


DEPTH IN METERS AT BOUNCE. INITIAL ELEVATION ANGLE = 30.00 


HEADINGS : 10.00 
BOUNCE 

3160.6 
4042.9 
4343.8 
3627.9 
2780.2 
2190.1 
1823.4 
1606.3 
1490.1 
1452.0 
1481.7 
1586.9 
1785.7 
2126.5 
2663.6 
3451.5 
4246.8 
(4261.3 
3486.1 
2713.1 
2187.8 
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TABLE VI 


RAY FAN SYMBOLS - used to denote location of bottom reflections in the 
ECTRACE plots. Thirteen different symbols are used before repetition. 
All rays of the same ray fan (initial elevation angle) will indicate a 
bottom reflection with the same symbol. Each ECTRACE plot lists the 
symbols used for the horizontal and vertical plots. 


a NR RNR $5 a 


INITIAL ELEVATION ANGLES 
RAYFAN SYMBOL THETA i 


1 ul -10.000 
2 O ~9.000 
3 ‘A -8.000 
mn 4 ~7.000 
S . -6.000 ; 
6 © -5.600 . 
7 As -4.000 
8 x -3.000 
| g Zz -2.000 
. 10 Y -1.000 
= 11 x 0.000 
z le * 1.000 
} 13 x 2.000 
7 1u =] 3.000 
15 © u.000 
16 A 5.000 
| 17 o 6.000 
18 x 7.000 
19 © 8.000 
20 4 9.000 
21 x 10.000 i 
HEADINGS 90.000 +/- 0.000 | 


TABLE VII 
ECCOM printed output. The rays which are plotted are listed, and their 
initial heading, initial elevation, and the number of points (initial plus 
reflections) are indicated. In this run the plot produced is shown in 
Fig. 7. 
ECCOM 51 *** SYNTHETIC BASIN *** 
SYNTHETIC SEA BED 
COMPUTER GENERATED BASIN 
XST XDV YST YDV CNX CNY CNLAT CNLON 
0.0 10.000 0.0 10.000 75.000 75.000 0.0 0.0 
INTIAL INTIAL 
ELEVATION HEADING POINTS SCALE FACTOR= 0.340 
40.000 40.000 21 
40.000 44.000 21 
40.000 48.000 21 
38.000 42.000 21 
38.000 46.000 21 
36.000 40.000 21 
36.000 44.000 21 
36.000 48.000 21 
34.000 42.000 21 
34.000 46.000 21 
60.000 40.000 21 
60.000 44.000 21 
60.000 48.000 21 
51.333 42.000 21 
51.333 46.000 21 
42.667 40.000 21 
42.667 44.000 21 
42.667 48.000 21 
34.000 42.000 21 
34.000 46.000 21 
60.000 5.000 21 
60.000 9.000 21 
60.000 13.000 21 
62.000 7.000 21 
62.000 11.000 21 
64.000 5.000 21 
64.000 9.000 21 
64.000 13.000 21 
66.000 7.000 21 
66.000 11.000 21 
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APPENDIX H 
FIGURES 

Program and Data Utilization Flowchart 

Contour Plot of Lofoten Basin (GENBOT) 
Horizontal Ray Plot Over Mohns Ridge (ECTRACE) 
Contour Plot of Mohns Ridge (GRDCHK) 

Complete ECTRACE Plot for a Synthetic Basin 
Contour Plot of a Synthetic Basin (GRDCHK) 

ECCOM Plot for a Synthetic Basin 
NRL N. Atlantic Contour Chart 
Perspective Surface Plot - Seamount (SYNGEN) 
Perspective Surface Plot - Basin (SYNGEN) 
Perspective Surface Plot - Mohns Ridge (SYNGEN) 
Contour Plot of a Synthetic Seamount (GRDCHK) 
Complete ECTRACE plot for a Synthetic Seamount. 
Horizontal Ray Plot over a Synthetic Seamount (ECTRACE) 
Vertical Ray Plot in a Basin (ECTRACE) 

Vertical Ray Plot over a Seamount (ECTRACE) 
Horizontal Ray Plot in a Basin (ECTRACE) 


ECCOM Plot for a Synthetic Wedge 
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Program and Data Utilization Flowchart. 


FIGURE 1. 
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Vixen ELEVATION ANGLES 
RATER OST MEOL THET! 


t 2 25.000 
2 s 30.000 
E a 8.000 


MEROINGS 50.000 °/- 30.000 


90 27.ca 5 x cry ans 08 CY) 90 
ROMNS AIOGE CENTER LAT/LON 71~37N/000-02" 


FIGURE 3. ECTRACE TRIAL 82C performed over a depth-grid matrix of 
the Mohns Ridge. In this run three ray fans are used 
with eight rays per ray fan. This may be overlaid on 
a contour plot of this area (Fig. 2). 
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FIGURE 4. A GRDCHK plot made on a matrix generate 
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ECTRACE TRIAL 67, performed in a synthetic parabolic basin 
(FIGURE 6). In this run eight ray fans are used with initial 
elevation angles varying between -30 and +30 degrees. Since 
there is only one ray in each ray fan, the initial heading 
for the rays will be 10 degrees. The vertical plot reveals 
CZ, sofar ducted, and SRBR rays. 


CCOGN $1 mee SCTNTHETIC G06lN | um 


| ¢ aa ae 
| . SYNTHETIC SBA SED CORPUTER GENERATED 


FIGURE 7. Convergence in a parabolic basin (Fig. 10). This is an 
ECCOM product using punched output from three separate 
ECTRACE runs. TABLE VII identifies the ray headings, 
elevation angles, and number of points. 
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FIGURE 8. 


A reduced representation of a chart titled "Bathymetry of 
the Norwegian~Greenland and Western Barents Seas" published 
by NRL (Washington, D.C.). The actual chart measures 4 by 
3 feet depicting finely grained contours. GENBOT uses 
bathymetry data files of this area to generate the matrix 
grids used by ECTRACE. 
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FIGURE 9. A perspective surface plot produced by G3DP of a conical 
seamount generated by SYNGEN. 
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FIGURE 10. A perspective surface plot of a synthetic parabolic basin 
produced by SYNGEN. 
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FIGURE 11. A perspective surface plot produced by G3DP of a depth 


matrix representing a portion of the Mohns Ridge, see 
Fig. 2. 


1.0 


frases at + edigshrale- 


ae! 


38,5 


113. 


151.6 


plbineent eden ee + 


: SONPUSEAS GENERAY ED SEAMOUNT 
FIGURE 12. A graphical analysis performed by GRDCHK of a depth matrix 
representing a conical seamount. The upper portion is a 
series of displaced plots of selected columns of the matrix. 
The lower portion is a contour plot. The seamount was 
created by SYNGEN, see Fig. 9. 
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An ECTRACE plot of two ray fans in a synthetic parabolic 


basin (FIGURE 10) generated by SYNGEN. 
initial elevation angle of 30 and 45 degrees (positive 


angle downward). 
initial headings 360 and 020 degrees, respectively. The 


plot was terminated after each ray had accumulated 23 


bottom bounces, 
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FIGURE 13. 
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FIGURE 17. An ECTRACE horizontal plot product of a single wide ray fan 
in a synthetic basin, All rays share the same bottom con- 
tact symbol. The number at the last symbol identifies the 
ray's order in the printed and punched history. 
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A horizontal plot produced by ECCOM of four ray fans propa- 
gating over a six degree wedge. The punched output of three 
ECTRACE runs was combined and then every other ray data his- 
tory card sequence was removed. There are a total of 40 rays 
depicted varying in initial elevation angles from 20 to 26 
degrees in two degree increments. The initial headings vary 
from 31 to 49 in two degree increments. The main purpose of 
ECCOM is to display the horizontal curvature effect of a 
large number of rays. 
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APPENDIX 0 
FORTCLGW JCL 


’ The job control language (JCL) card for FORTCLGW used in several of 
the programs was 


//EXEC FORTCLGW.REGION.GO#250K . 
This JCL card identifies and calls the following JCL: 


XXFORTCLGW PROC DEST=A, IMSL=SF 

XXFORT EXEC PGM=IEYFORT,REGION=180K 

XXSYSPRINT DD SYSOUT=&DEST , DCB=(RECFM=FBA, LRECL=120 , BLKSIZE=3360) 
SSSYSLIN DD UNIT=SYSDA,SPACE-(CYL, (1,1)) ,DISP=(,PASS) ,DSN=&&SYSLIN, 
XX DOE=(RECFM= FB, LRECL=80 ,BKJSUZE=800) 

XXLINK EXEC PGM=IEWL,REGION=180K,PARM='MAP,LIST' ,COND=(4,LT) 
XXSYSLIB DD DISP=SHR,DSN=SYS1.FORTLIB 


f XX DD DISP#=SHR,DSN=SYS1.MPSLIB 
XX DD DISP=SHR,DSN#SYS3.IMSL.&IMSL 
’ XX DD DISP=SHR,DSN=SYS1.VTECPLOT 
XXSYSLMOD DD DSN=&T (PROGRAM) ,UNIT=SYSDA, DISP=(,PASS) , 
XX SPACE=#(CYL,(3,1,1)) . 
XXSYSLIN DD DSN#*,FORT.SYSLIN,DISP=(OLD, DELETE) 
XX .D DISP=SHR ,DSN=SYS1.PROCLIE(VMAPP) 
XX DD DDNAME=SYSIN 


XXSYSPRINT DD SYSOUT#&DEST,DOE=(RECFM=FBA, LRECL=121,BLKSIZE=1210), 
XX SPACE=(CYCL, (1,1)) ,UNIT=(SYSOUT ,SEP=SYSLIN) 
XXSYSUT1 DD UNIT#=(SYSDA,SEP=(SYSLIN,SYSLMOD,SYSLIB)), 
XX SPACE#=(TRK, (10,5)) 

XXPLOT EXEC PGM=IEVMAPP ,COND=(4LT) 

XXSTEPLIB DD DSN#=SYS1.VTECPLOT ,DISP=SHR 

XXPLOTLOG DD SYSOUT#=A 

XXVECTR1L DD DISP=(OLD,DELETE) ,DSN=&&VECTR1 
XXSYSVECTR DD DISP#=(OLD,DELETE) ,DSN#&&VECTR2 
XXSYSVECTOR DD SYSOUT#=(A, ,5555) 

XXVECTTAPE DD DUMMY 
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